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FOREWORD 


The report discusses the development of a systematic 
methodology to aid pilots in accurately following aircraft 
flight test trajectories. Modern control methods are used to 
compute time histories for displays which both reduce pilot 
workload and make test trajectories consistently accurate. 
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SECTION 1 


INTRODUCTION 


It is often necessary to follow desired time histories of 
several variables during an aircraft flight test. The pilots 
often find it difficult to follow such trajectories accurately, 
leading to increased test time and pilot workload. This report 
shows that modern control techniques can be applied to signifi- 
cantly improve the accuracy with which complex trajectories can 
be flown. In addition, the pilot inputs are concentrated below 
2 rad s - ^ reducing pilot workload and fatigue. 

The study involves computing three display variables on 
which pilot pitch, roll and throttle commands can be based. The 
pilot perceives a three axis uncoupled control problem because 
all couplings are handled by the flight trajectory controller. 

The emphasis has been placed on the development of the 
techniques. An F-15 nonlinear simulation has been used for 
demonstration. The techniques are general and should be appli- 
cable to various combinations of trajectories, aircraft and 
displays . 


3 . 1 SUMMARY OF APPROACH 

The principal steps for systematic development of flight 
trajectory controllers are 

1. Development of the flight test plan, 

2. Control design model specification, 

3. Controller design, and 

3. Controller validation by analytical simulation 
and pi 1 ot-in-the-loop-simulat ion . 



The flight test plan prescribes the aircraft with all subsystems 
configured for the given test, the trajectory or classes of tra- 
jectories to be flown, and the pilot's tasks during the specified 
trajectory maneuvers. Figure 1-1 depicts the design procedure 
after the flight maneuver and system models have been specified. 

1.2 RESULTS 

The major result from this study is a set of tools suitable 
for developing flight test trajectory controllers. Attention has 
been given to keep the techniques as general as possible. The 
techniques should have a general application to a wide class of 
problems where quantities displayed to a pilot should be appro- 
priately computed to improve pilot effectiveness and to reduce 
workload and fatigue. 

The development of the control laws is straightforward 
using an advanced control design program once the model and the 
performance index have been specified. Thus, major part of the 
attention must be given to these two tasks. Special linearization 
tools were developed for the F-15 simulation to develop a suitable 
model for control. Integration of pilot model is a difficult task 
because the pilot model itself depends upon the trajectory con- 
troller. This was achieved by prespecifying a suitable represen- 
tation for the aircraft and controller dynamics and basing a 
pilot model on this dynamics. It is necessary to roll-off the 
control law at high frequency, and frequency shaping methods are 
essential to achieve necessary gain margin at high frequency. 

Pilot-in-the-loop and flight test validation of controllers 
are needed. It may be desirable to integrate pilot inputs with 
autopilots to achieve even better overall performance. 

1.3 REPORT ORGANIZATION 

This report is organized as follows: Sections 2 through 4 

describe systematic methods for designing flight test trajectory 
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Figure 1-1 


Flight Test Trajectory Control Design Procedure 
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controllers, illustrated on the F-15 aircraft. Section 2 cate- 
gorizes how the equations are developed for different types of 
maneuvers. Section 3 develops linear models for each major 
subsystem, and Section 4 describes the control design methods. 

Section 5 describes the control law synthesis for the zoom 
and pushover maneuver, and linear and nonlinear evaluation of 
the control law. Section 6 documents this control algorithm for 
use in a batch simulation, manned simulator or flight test. 

General conclusions are given in Section 7. 


i 
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SECTION 2 


TRAJECTORIES FOR FLIGHT TEST MANEUVERS 


A large class of trajectories can benefit from flight test 
trajectory controller design. The model for the control law 
and the development of the control law will depend strongly on 
the particular trajectory to be flown. 

Aircraft flight test trajectories could be based on inertial 
reference (e.g., level-turn or 3-D guidance) or reference with 
respect to another vehicle or vehicles (e.g. in air-to-air combat). 
To place the flight test trajectory control design problem in 
a proper framework, both of these cases will be discussed. 

2.1 SINGLE VEHICLE FLIGHT TEST TRAJECTORY 

We divide single vehicle flight paths into those which 
require continuous control along the trajectory and those that 
specify a final flight condition. In either case the flight 
test trajectory could be specified in terms of one of the 
fol lowing : 

1. constraints on position components, 

2. constraints on velocity components and altitude, 

3. constraints on combinations of load, speed and 
altitude . 

Combinations of these constraints could also be considered. 

Constraints on Position Components 

Examples of trajectories which involve position constraints 
along the flight path are 
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1. 4-D guidance (x(t), y(t), h(t) are given functions of 
of time) . 

2. 3-D guidance (x, y, h are related to each other, 
e.g., fly along a hypothetical wire in space). 

Examples of 3-D guidance are approach to landing, 
terminal area flight paths and threat evasion for 
reconnaissance aircraft and bombers. This also 
includes straight and level flight and flights 
along predetermined paths. 

Examples of trajectories which specify position constraints at 
final trajectory point are: 

1. 4-D specification (arrive at a certain point, at a 
certain time, e.g., touchdown on the runway at a 
specified point). 

2. 3-D specification (fly-to-VOR, terrain following). 

Note that each of these trajectories requires position measurement. 
The 4-D guidance trajectory indirectly specifies velocity and 
acceleration components. Thus specification of position compon- 
ents is the most comprehensive constraint on the trajectory. 

Such a rigid constraint may be unnecessary for most test maneuvers. 

Constraints on Velocity Components and Altitude 

While the horizontal position components do not, in general, 
affect aerodynamic variables, the altitude determines density 
and by itself affects dynamic pressure. Thus it must always be 
considered as a possible variable to be constrained. In fact, 
the altitude and dynamic pressure are so important that a 
majority of flight test trajectories will define the altitude 
profile (this includes maintaining constant altitude). 

Examples of this class of trajectories are: 

1. u(t), v(t), w(t) and h(t) [in other words, Mach 

number, dynamic pressure, a(t) and 8(t). 6 ( t ) may 

be zero. 

2. Mach number, angle-of-attack and dynamic pressure 
(as in shuttle tile tests). 

Various other combinations of velocity components and altitudes 
could also be specified. 


6 



Mach number, angle-of -att ack and altitude constraints could 
also be desired at one point on the trajectory. Zoom and pushover 
is a trajectory where angle-of-att ack , Mach number and altitude 
are specified at one point on the trajectory. Note that since 
sideslip angle is assumed zero, zoom and pushover requires a 
certain set of u, v and w at a certain altitude. 

Constraints on Combinations of Load, Speed and Altitude 

The trajectory specifications could involve components of 
loads along the three axes, velocity components and altitude. 

The typical load specification will consist of desired vertical 
acceleration. The desired value of the lateral acceleration 
is usually zero. The total speed is specified in lieu of the 
fore-and-aft acceleration. 

Many combinations of load, speed and altitude specifications 
are possible. Some examples are as follows: 

1. A constant load, constant Mach number level turn, 

2. A constant Mach number, constant altitude spiral turn. 

2.2 MULTIPLE VEHICLE FLIGHT TEST TRAJECTORIES 

Often the desired flight trajectory for an aircraft depends 
upon the position and flight test trajectory of other vehicles. 
Typical examples are collision avoidance or air combat with air- 
to-air missiles. The specification is typically based on the 
position of a target aircraft with respect to the aircraft whose 
trajectory is being controlled. 
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SECTION 3 


SYSTEM MODEL FOR TRAJECTORY CONTROLLER DESIGN 


The overall system model for the flight test trajectory 
control analysis has the following components: 

1. aerodynamics and kinematics, 

2. hydro-mechanical control system (MCS) and command 
augmentation system (CAS) including actuator models, 

3. engines, 

4. sensors, and 

5 . pilot . 

The flowchart of the complete model is shown in Figure 1. 


Figure 3-1 Schematic of System Model for Flight Test 
Trajectory Control Analysis 
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3.1 AERODYNAMICS AND KINEMATICS 


The structure of mathematical equations governing the motion 
of a rigid aircraft may be derived from kinematics. Aerodynamics 
will be linearized, as necessary. A different approach must be 
used to smooth out any discontinuities in slope. 


Consider an aircraft state vector 
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(3.1) 
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and 


(3.2) 


The equations of motion can be written as 


X = f(X,u) 


(3.3) 


A detailed nonlinear aerodynamic model of the F-15 from 
the NASA Dryden Flight Research Center Simulation SIMII was 
linearized by 

a) trimming the aircraft at the desired flight 
condit ion , 

b) perturbing the six aerodynamic force and moment 
equations with respect to all dynamic variables of 
interest, and 

c) projecting these stability derivatives onto both 
body and stability axes (see Tables 5-1 and 5-2 in 
Subsection 5.2). 
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3.1.1 Longitudinal Aerodynamic Model 


The simplified linear longitudinal equations developed 
further in Section 5.2 for the zoom and pushover maneuver are: 
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- g //sin a Q - X a + g 
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where 


6V - speed 

6a - angle-of-attack 
q - pitch rate 
69 - pitch attitude 
6h - altitude 


6g B - stabilator deflection 
5T/m - thrust acceleration 


are perturbation states about the trimmed condition in stability 

axes . 


3.1.2 Coupled Lateral /Longitudinal Equations 

For high bank level turns trimmed at a constant turn rate, 


p o = - K sin e o 

q o = / cos 0 O sin <P Q 
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r o = *o cos 0 o cos *o • 

The coupled linear equations in body axes are 
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where r x = (I y -I z )/I x , Q y = <I Z -I x >/I y > P z = 

3.2 MCS AND CAS MODELS 

The mechanical control system (MCS) and control augmentation 
system (CAS) use pilot commands and instrumentation system measure- 
ments as inputs. The outputs of this system are the aerodynamic 
control surface deflections and the thrust command. 

Only those dynamics in the pilot's range of attention and 
response need to be modeled (.1 to 20 radians per second). 

Symmetric simple nonlinearities such as position and rate limits 
have no memory and should be avoided in trajectory controllers. 

An overview schematic of the longitudinal MCS/CAS system 
is given in Figure 3-2. The mach servo loop is fixed at a con- 
stant gain obtained from SIMII at the trimmed condition. Des- 
cribing function analysis indicated that the nonlinear control 
elements could be realized at the zoom and pushover trimmed 
condition as shown in Figure 3-3 (see Section 5.2 for the linear 
equations) . 
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Figure 3-3 


F-15 MCS/CAS Longitudinal Linear Model (Roll, Pitch, Yaw CAS on) 









3.3 ENGINE MODEL 


Simple as well as highly sophisticated models have been 
developed for engines such as F100-PW-100. A simple first-order 
lag model should suffice for the purpose of the study. The engine 
thrust T is related to the commanded engine thrust T^, as 

T = _ I t + - T_ , (3.4) 

T T C 

where x is the time constant in the first-order lag model. 

3 . 4 SENSOR MODEL 

Typical instruments on board a test aircraft include a three- 
axis accelerometer, a three-axis rate gyro, angle gyros, angle-of- 
attack and sideslip angle vanes, altitude, rate-of-climb and total 
speed . 

We could use a model with first-order lag and random noise. 
Thus, the measured quantity y m will be related to the true y 
as follows 

: v = y - v + n (3.5) 

■ m ' m 

where t is the time constant and n is the random noise. All 
sensors give sampled measurements (with an assumed rate of 10 Hz). 

3.5 PILOT TASK DEFINITION 

The trajectory models described in the previous section 
specify a portion of the controller task; however, since the 
pilot is actually part of the overall controller his precise 
role, the variables that are displayed to him, and the controls 
he uses to zero the display must be considered. Alternative 
roles for the pilot in the overall controller include: 1) using 

the pilot only in the loop, nulling the displays driven by the 
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trajectory controller, 2) giving the pilot vernier controls to 
refine the performance of the autopilot, or 3) flying the 
maneuver automatically after the pilot has achieved the appro- 
priate inital conditions. 

The first role is desirable since the pilot's isolation of 
faults and corrective action is much more reliable than that of 
the autopilot. Only this first role, with the pilot fully in 
the feedback control loop will be used in this study. 

Besides the pilot role, the control design must also address 
which error signals should be displayed and which controls are 
most effective for a given trajectory. To simplify the pilot 
task the errors displayed should be dynamically decoupled with 
respect to the controls as much as possible, and also separated 
in frequency where possible. This enables the pilot to more 
easily coordinate and sequentially correct the trajectory errors. 

There are two other aspects of the pilot task that effect 
the controller design, the pilot's preference for the "plant 
controller" dynamics, which is the combined plant that he con- 
trols, and the pilot's preference for the speed of the closed- 
loop system; i.e., the speed at which he prefers to control. 

We want to optimize the pilot performance by giving him a system 
which he can best control. Since the pilot responds differently 
to different systems, iterative controller design could be used 
to converge on the combined pilot controller-plant model. 

For this study, a published model of the human operator controlling 

o 

a 1/s plant [1] will be used. Further research is currently being 

conducted, employing improved representations of pilot behavior 

using frequency shaping methods. Since the human operator adapts 

2 

to the given system the 1/s model will be a conservative des- 
cription of the pilot's actual performance in flight tests. 
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Pilot Model 


A pilot model presented in previous studies has been 


PCs) 


EQ (s 


e 


- TS 



+ 1 


(3 .6) 


where P^Cs) is the human equalization network, t is a 
neurological motor delay, and x an observation delay. Figure 3-4 
shows what has been called the control theoretic model of optimal 
human behavior. 


Pilot Modeling Difficulties 

The above control theoretic model requires the propagation 
of delay-differential equations which introduce considerable 
compl icat ion . 

Between the two delays and the control and filter weights 

in the optimal regulator and estimator, there is considerable 

flexibility in this model to fit empirical data. If the system 

is driven with sinusoidal inputs then the assumption that the 

human can perfectly compensate for a delay will be self-fulfilling 

as the human operator settles into a steady-state after observing 

several cycles. Time domain approaches which include random 

inputs with a band of frequency content would be more desirable 

and should indicate a much lower gain at high frequencies, rather 

than primarily the large phase lag. Finally, the model is valid 

2 

for a particular plant, in our case 1/s , when the dynamics are 
actually higher order and not a pure integrator in the mid-fre- 
quency range where the pilot responds well. 

The significant aspect of the human pilot, however, is his 

ability to identify and adapt to the plant dynamics, therefore 

2 

our assumption of a 1/s plant is a conservative one for the 

overall trajectory controller and pilot system. Thus, the given 
2 

1/s model will suit our purpose in developing trajectory con- 
troller design techniques. The frequency response shown in Ref. 
fl] is given in Figure 3-5. 
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The pilot limitation parameters are given in the following 
table: 


Motor Delay 

— 

t n 

100 msec 

Observation Delay 

T 

210 msec 


The pilot transfer function is 


P EQ (s) (9.52 - s) 
(s + 10) (s + 9.52) 


(3.7) 


A reasonable linear approximation to the above frequency response 
i s 
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3.6 SUMMARY 


A state space linear perturbation model for the overall 
system shown in Figure 3-1 appears as 


(3.10) 




where the states, controls, and noise inputs are: 


x L = location coordinates (x, y, z) 
x A = aerodynamic /kinematic states 
X MC ~ MCS/CAS states 
x £ = engine states 
x p = pilot states 
x = sensor states, 

b 




vertical display error 
horizontal display error 
Mach display error 


w = gust inputs 


n s = sensor noise 
n p = pilot motor noise 
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SECTION 4 


CONTROL DESIGN PROCEDURE 


The overall control-design procedure consists of three major 
steps : 

1. Model reduction : Development of control-design 

models from detailed models of the previous section. 

2. Control design : Design of a controller/f ilter or 

output-feedback logic to specify the control 
structure . 

3. Controller simplification : Reduction of the con- 

troller obtained in the previous step to simplify 
implementation and for more robustness. 

The steps in the overall procedure must be carefully inte- 
grated to support assumptions and approximations made in each 
of the steps. The development of specific design methodologies 
which support this integration is the basis of our control design 
methodology for flight test trajectory controllers. 

We discuss each of these steps in some detail . 

4.1 MODEL REDUCTION 

Model reduction is required for robustness and to simplify 
control law computation and implementation. 

Criteria and computation procedures for model reduction 
are discussed in this section. The reduced model must contain 
essential elements of the dynamics such that the resulting 
control law meets the design requirements. Roles of poles, zeros 
and residues of transfer functions in closed-loop control design 
are discussed first followed by selection of specific criteria 
and numerical procedures. 21 



Poles, Zeros and Residues : In linear systems, the closed-loop 

behavior is dictated significantly by several transfer function 
properties. Three of the most significant are the poles, the 
zeros, and the residues. 

Consider a linear model with nxl state vector x, qxl control 
vector u, and pxl output vector y: 

x = Fx + Gu , 
y = Hx + Du . 

The transfer function between y and u is 

y(s) = [H(sI-F) _1 G + D] u(s) 

A T(s) u ( s ) . (4.3) 

The roles of transfer function properties are best explained 
by considering a single-input single-output (SISO) system 

y(s) = T ( s ) u(s) (4.4) 


(4.1) 

(4.2) 


The characteristic values , i = 1 , 2 , . . . , n of the 
T(s) denominator are the open loop system eigenvalues and 
indicate system stability properties. The transfer function 
may be written in terms of zeros z^, i = 1 , 2 , ..., n or 

residues r^ , i = 1 , 2 , ..., n as follows: 


K(s-z 1 ) (s-z 2 ) ... (s-z m ) 

T(s) = — 7 r — r -7 - — r r — r- : Zeros Representation (4.5) 

(s-XjMs-lg) ••• (s-A ) 


i=l 


R 


-*i 


: Residue Representation (4.6) 
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Consider now a feedback matrix C(s) with gain a . The 
closed-loop transfer function is 


T c (s) 


T(s) 

1 + aC(s) T(s) 


( 4 . 7 ) 


which may be written as 


T c (s) = 


K(s-z- ) ... (s-z— ) 

1 m 

(s-X ]L ) (s-A 2 ) ... ( s — A— ) 


It is easy to see that for smaller a , 


dX i 


3 a 


cx ->-0 


= cxC ( X i ) r ± 


( 4 . 8 ) 


( 4 . 9 ) 


Thus r^ dictates the behavior of the pole for small gain. When 
a is large, the finite closed-loop poles are the zeros of C(s) 
and T ( s ) . 

To summarize, the residues of the transfer function describe 
the low-gain properties and the zeros the high-gain properties. 
Both zeros and residues are important in closed-loop control 
design. Therefore, an ideal reduced-order model should maintain 
the residues of the retained poles and zeros in the spectrum of 
interest. Unfortunately, both zeros and residues cannot be 
preserved simultaneously. The attempt in reduced-order modeling 
methods is to maintain either poles or zeros or provide approxi- 
mations to both of them. 


Retention of Residues 


The procedure for retaining residues is implemented as 
follows. Let F, G, H and D be in modal form and assume 
that the first n modes are retained. The state equations are 
then written as 
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+ Du . 


(4.11) 


y - («! | h 2 ) 

Let the average frequency of the retained model be w . The 

Bi 

state Xg is approximated by 

x 2 = Re[( ju a I-F 22 ) _1 G 2 ] u , (4.12) 

where Re(«) represents the real part of (•). Therefore, the 
reduced model is approximated by 

x^ = F^x^ + Gu , 

y = H 1 X 1 + tD +Re[H 2 (:ia) a I “ F 22 ) ” lG 2 ]] u * (4.13) 

I f we are interested in matching the dc gain of the reduced model 
to that of the high-order model, oj may be set to zero. Note 

BL 

that only the real part of the term in Eq. 4.12 is retained 
because we want reduced-order models with real coefficients. 

Retention of Zeros 



Suppose we want to retain the first m zeros and n poles. 
Let ~ be the average frequency. The transfer function of 

BL 

Eq . 4.5 is simplified to 



(s) 


g(ju a- 2 m+l } ( jy z M } < s ~ z ! ) ••• (s ~ z m ) 
( j u 'a" X n+l ) *•* ^ ••• ^ S-X n ^ 


(4.14) 


The first term is approximated by a real gain. Each input -output 
transfer function is simplified as in Eq . 4.14. These simplified 
transfer functions are reconverted into a st ate-variable 
description . 
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The computation time requirements to obtain reduced-order 
models which retain zeros are very high. Poles and residues 
of high-order models are determined much more easily than zeros. 
This model-reduction procedure, therefore, has been used rarely. 
Methods that maintain system transmission zeros have been des- 
cribed using zero directions [2]; however, a systematic design 
method employing MIMO zeros has not yet been developed. The 
transmission zeros are the frequencies at which there is no 
output to any combination of inputs. They are the eigenvalues 
of the generalized eigenvalue problem 


X 

= z 

' I 

O' 


X 

u 


0 

0 


u 


(4.15) 


Mode Detection Criteria 


If the model is to be reduced to n' states, the number 
of poles will drop to n'. Criteria have to be developed to 
determine the set of poles which may be dropped, along with the 
corresponding set of modes. These criteria must consider the 
fol lowing: 

1. Any mode which is uncontrollable and undisturbab le 
or unobservable can be dropped. Either of these 
conditions corresponds to a zero residue or a 
perfect pole-zero cancellation in all transfer 
funct ions . 

2. Controllable modes should not necessarily be dis- 
carded even if they are nondisturbable and add 
nothing to the cost functional, since they may be 
excited by the control actuators. 

3. Highly controllable modes in the bandwidth of 
interest should be retained even if they are not 
disturbable or observable. Similarly, highly 
observable modes should be retained for robustness 
reasons . 

4. Proper mode ordering based on cost analysis [3] 
should include the above factors in addition to 
performance considerations. 
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It is clear that many simplifications of the model form 
are possible. In general, the analyst would start from the 
simplest model and carry the procedure through to evaluation. 

If the performance is unsatisfactory, a more detailed model 
should be considered. 

4.2 SELECTION OF CONTROL LAW AND STATE ESTIMATOR 

Multivariable linear-quadratic-regulator (LQR) control 
theory is a powerful tool for the development of flight test 
trajectory controllers. If the measurements have unacceptable 
noise or if unmeasured states are needed for feedback, a state 
estimator may also be designed to supplement the regulator. 

4.2.1 Basic Linear Quadratic Regulator Approach 

The basic approach is described first. Extensions of the 
approach for the test trajectory control problem will then be 
covered. Consider a linear system model in which the states 
x and inputs u represent perturbation around nominal values 
x Q and u Q 

x = Fx + Gu . (4.16) 


In this problem, u represents the variables to be displayed 
(nominally roll, pitch and Mach number errors) and x depends 
on the particular model dynamics selected (u, w, z, T and 
h m for the simplest model). The control design problem selects 
u to drive the perturbation states x to zero. The LQR formu- 
lation optimizes a quadratic performance index of the form 


J 



(x^Ax + u^Bu) dt 


(4.17) 


The control law is obtained by solving the Riccati equation 


- A + SGB 1 G T S 


26 


S 


SF - F T S 



and 


u ( t ) = - B VS x ( t ) = C(t) x ( t ) . (4.18) 

To simplify implementation, a steady-state solution could be 
used where the S is set to zero and a constant value of the 
gain matrix C(t) is computed. This constant control gain can 
be used throughout the maneuver. Note that since x and u 
are perturbation values, the steady state x Q and u Q must 
be subtracted from the measured quantities prior to computing 
the control input, and is added into the final control applied. 

Each maneuver is specified in terms of desired values of 
the 12 aircraft state variables or functions of these variables 
(the 12 variables are the three components each of position, 
velocity, Euler angles and angular rates.) The specification 
can be in one of the following forms for each state (or function 
of states): 

1. desired constant value, 

2. desired constant value for integral or derivative, 

3. desired time history, or 

4. desired value at some point on the trajectory. 

4.2.2 Constant Values 


A desired constant value is equivalent to having trim values 
as discussed above. The steady-state states and controls neces- 
sary to achieve the desired outputs can be solved from the equation 



(4. 19) 


If the system has open-loop zeros at the origin, this inverse 
does not exist, or in other words, there is no neighboring 
steady-state value. In this case integral error states can be 
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added to achieve the desired constant value. Constant integral 
as wel] as derivative states are discussed next. 

4.2.3 Constant Integral or Derivative Regulation 

Integral control ; To account for errors in computation 
of trim conditions and other low-frequency modeling errors, it 
is necessary to place a penalty on the integrals of the error 
signals. For example, if a state z is being held at z Q , we 
can define 

i = z - z o (4.20) 

and place an additional penalty on £ as follows: 

C T . (4.21) 

Regulation of derivatives of states : To control derivatives 

of states, additional terms are added to the performance index. 
For example, to make Ex follow E Q , the following term is 
added to the integrand of J : 

(Ei - E q ) T A e (Ex - E q ) , 
or 

[ E ( Fx + Gu) - E o ] T A g [E(Fx + Gu) - E Q ] . (4.22) 

The control laws can be computed directly for the new performance 
index and will be of the form: 

u ( t ) = - Cx(t ) + C e E Q . (4.23) 

(References [4] and [5] describe a new stable and reliable algo- 
rithm for solving the algebraic Riccati equation with cross- 
coupling between state and control cost.) 
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Integral control is valuable for zero steady-state errors 
although it can lead to a sluggish controller if the feedback 
consists primarily of the integral state, since it takes time 
for the error to build up. 

Derivative control on the other hand gives excellent response, 
and in fact, is called dynamics matching in the context of implicit 
model following. A typical trajectory where such a control law 
would be possible is the specification of loads, where linear 
combinations of the derivatives of rate equations are to be 
regulated about a constant value. 


4.2.4 Controlling the Time History of a State Variable 

The technique can be extended to make certain variables 
(say Tx) follow a time history T(t). The integrand of the 
performance index is extended as follows: 

(Tx - T ( t ) ) T A (Tx - T (t ) ) . (4.24) 


The control law now requires the solution of a differential 
equation backward in time. 

The complete problem formulation as described in Ref. [6] is 


where 

x 


1 

2 



- H o x(t o 




[ Z 



+ 


U 


2 

_l J dt 
B 



H f x(t f ) 



(4.25) 


Fx + Gu 


(4.26) 


given Zq , H Q , A Q , z(t), H(t), A(t), B(t), z f , H f , A f , 

F(t), G ( t ) , t Q , and t f . 
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x(t) is the state vector of a dynamic system. 

H(t) x(t) is a vector output history. 

z(t) is a desired vector output history (prespecified). 
u(t) is a vector control input history. 

H q x ( t q ) and H^x(t^) are initial and final vector outputs. 
z Q and z f are desired initial and final vector outputs. 

Aq , A f , A(t), and B(t) are relative weighting matrices 

to be specified by the control designer. 

The desirable solution for this "follower" problem is a backward 
sweep solution derived from the Euler-Lagrange equations, of 
the form 


A ( t ) = - }, B (t) + S B (t) x(t ) . (4.27) 

Differentiating (4.27) with respect to time and using the Euler 
Lagrange equations yields the backward sweep equations: 

- s B = s b f + f t s b - s b g“ 1 g t s b + H T AH , 

S B (t f ) = HjA f H f , (4.28) 

- \ B = (F - GC b ) T A b + H T Az , 

^ B (t f ) = H f A f z f , (4.29) 

where 

C B = B _1 G T S b . (4.30) 

Equation (4.28) is a Riccati equation for the symmetric matrix 
Sg. It must be integrated backward from t = t^ to t = tg 
along with (4.29) for the vector Ag. 

-1 T 

If Cg(t) and B G Ag(t) are stored on the backward sweep, 
then x and A can be integrated forward to determine u(t): 
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X 


Fx + Gu , 


(4.31) 


where 


u = U B - C B x , (4.32) 

u B = B _1 G T X b . (4.33) 

Note that this set of equations involves linear feedback of the 
state vector x(t). Equation (4.32) is a feedback plus feed- 
forward from of solution. 

For linear time-invariant systems with which we have approxi- 
mated our system, the Riccati equation need not be integrated 
but rather a vector forcing function using model decomposition 
for the linear two point boundary value problem (TPBVP) (see [7] 
for the distinct eigenvalue case and [8] for the general eigen- 
structure case). However, except for 4-D guidance, we do not 
need to control with a finite time formulation, but rather can 
employ steady-state control laws. This is equivalent to saying 
that z(t f ) persists as a constant value for a long period. 

This further simplification gives a constant gain solution 
for the control 


u = B 1 G T A b - B 1 G T SgX = u B - CgX , and (4.34) 


Sg(t) = Sg = AX 


-1 


where 


X. 

A 


(4.35) 


is an orthogonal basis for the stable eigenvalue of the Euler- 
Lagrange equations derived from the Hamiltonian. 


Three approaches to solving these equations without inte- 
grating Ag(t) backward in time are: 

1. Take Ag ^ 0 which implies that Ag % - F £ H Az , 
so that u = Ug - CgX, can be computed from 
z and x . 
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2. Consider z(t) as a constant, so that 
u = u Q - CgCx - x Q ) . 

3. Approximate z(t) as a function of the states 

a. A linearization of a nonlinear relationship 
between the states with integral error to 
correct the linearization. 

b. A random walk which can be integrated forward 
in time with x. 

1. and 2. above are equivalent, both relying on the assumption 

that the rate of change of z(t) is much less than the slowest 

T 

eigenvalue of F = (F - GCg) . This is a reasonable assumption 
for most trajectories, since commanded outputs are not rapidly 
changing and the controllable roots of F £ can be made faster 
with feedback. 

To approximate z(t) as a function of the states, the first 
alternative is not particularity attractive since an accurate 
linearization may not exist and the integral error tends to slow 

down the response and may also constrain the trajectory unneces- 

sarily. The random walk is a much more desirable way to generate 
a linear trajectory. The random walk also is affected by the 
rate of change of z(t) but not to the same extent as considering 

dz/dt =0 or lg = 0. The overall variance of the error in 

y may be comparable to the other approaches, but the model can 
be "tuned" to give a smaller spectral density in the range of 
i nterest . 

Changing the direction of integration in Equation (4.29) 

gives 

= (F - GCg ) T Ag + H T Az( t ) , (4.36) 

with an unknown initial condition. 

The appropriate initial condition Ag(t 0 ) appears as an 
impulsive input on the above equation. This is true at any 
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initial condition as well, hence we can merely append a white 
noise source to equation (4.36), an additive random walk term 


L = (F - GCtj) T A„ + H T Az(t) + n. , with 
n d n A 

X n (t Q ) = 0 , (4.37) 

n A = N(0 ’ V * 

and Q n can be chosen to bracket the desired value for Ag(t oo ) . 
Since 


>B (t «) - « Az(t.) , 


while Eq. (4.37) gives 




and Eq . (4.37) implies 


0 


'Xf m 

A + A F + Q 
n n c ^n 


=> 



(4.38) 


The RMS values from A„ give the interval about A (t) which 
should bracket the theoretical steady-state value Agtt^) with 
an acceptable level of confidence or integral error states can 
be added for A . 


More simply, the concept of the costate, A can be 

D 

dropped by using several appended random walk states driven by 
( Hx ( t ) - z(t)). 


4.2.5 Controlling a Set of States at Some Point on the Trajectory 

There are several ways to solve this problem. One approach 
is to select a state time history which leads to the desired 
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tra.iei tory point and then to derive a control law which follows 
the time history. It is also possible to define the trajectory 
point as a function of several states and then to hold this 
desired state fixed while other states are varied. The particular 
selection will depend on the desirable nature of the intermediate 
t ra ject ory . 

A combination of the above controllers can be used, if 
necessary . 


4.2.6 Extensions and Specializations 


Several extensions are proposed to improve the quality of 
the control laws. These extensions are designed to minimize 
the effects of modeling errors, lack of needed measurements and 


1. Estimator : If all states are not measured or if 

there is excessive random or quantization noise, 
v, a state estimator can be used. If y is a 
vector of measurements which is linear in the 
states , then 

y = Hx + v . (4.39) 

An estimate of the state can be obtained by 

x = Fx + Gu + K(y - Hx ) , (4.40) 

and K is obtained by solving a Riccati equation. 

(See the next section on controller simplification 
for a desirable algorithm in the case where some 
but not all of the measurements are noise free.) 

2. Frequency-shaped control law ; As shown in Appendix 
A, the control weightings can be made functions of 
frequency. By making the control weighting matrix 
an increasing function of frequency and the state 
weighting a decreasing function of frequency, the 
control activity can be concentrated in the mid- 
frequency range. The display will then have a 
minimum of high-frequency components. This shaping 
will also reduce the effect of neglected states 
with high natural frequency [9]. 


34 



4.3 CONTROLLER SIMPLIFICATIONS 


The controllers based on the LQR approach are simple and 
usually do not need to be simplified. LQG controllers involve 
a dynamic state estimator and are generally of the form 


Cx + D y . (4. 41) 

c c J 

The controller, therefore, has the same form as the system itself, 
except that: (1) the roles of input and output are switched, (2) 
the controller state vector depends on the choise of frequency 
shaping, and (3) K and C are known exactly (unlike G and 
H in the design model). Controller simplification is, therefore, 
similar to model reduction (see Section 4.1). 



Noiseless Measurements and Reduced Order Compensators 

The combination of the regulator and estimator can be viewed 
as a compensator for the open loop plant, where the compensator 
transfer functions are given by 


u ( s ) = T c (s) z ( s ) , 

u ( s ) = - C (si - F + FC + KH) -1 Kz(s) (4.42) 


This transfer function can be written in modal form as a parallel 
bank of first and second order equations. 


If some of the measurements are free of white noise, a 
r educed order compensator can be determined. Partitioning the 
measurements sc that the white-noise- free ones are at the bottom, 


T 

z = 


T 

z., 



the computed estimator gains are K = [K^ 01 . 
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These gains result from the reduced order Ricvati solution 


P = 


P 11 0 


0 


0 


(4.43) 


which can be solved for directly with recent stable and reliable 

a 1 gor i t hms [4,5] . 

The compensator feed-through term, , where 

u ( s ) = [ -C ( si - F + GC + KH) -1 K + D ] Z(x) , (4.44) 

can be computed from Zg which has no error, by considering 
that 



Z 1 


[ H ll 


Z 2 


*2 

The requested 

feedback 


C 

C D 

+ C 

with C D 

that portion 

z 2 = H 2 X " 

C D 

is the 


x + v 


and u = -Cx 


(4.45) 


the row space of C. If the intersection is zero there is no 
feed through, if the intersection is all of C then the dynamic 
compensation portion C v is zero. The feed-forward portion, 


u(s ) 

feedthrough 


= -C n H^z, = [ 0 *. D„ ] z , 

'2 


(4.46) 


is solved from the linear equations, 

,T 


h t n T 
H 2 D c r 


= C 


D 


(4.47) 
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All unobservable or uncontrollable modes of the controller 
should be dropped, since they have no effect on controller perfor- 
mance. In addition, poles with small residues can also be dropped. 
Note that the residues are physical matrices representing force 
or moment applied per unit error. 

Further selection of modes or states which should be retained 
in the low-order controller is a difficult problem. One approach 
is to drop one or two states at a time to determine the set of 
states which may be eliminated without loss of performance. To 
minimize computation time, a stepwise search procedure is needed. 

A modal cost function could be used in the search [10]. 
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SECTION 5 


ZOOM AND PUSHOVER ALGORITHM DEVELOPMENT 


This section describes the maneuver modeling, development 
of a reduced order linear design model, control design results, 
their evaluation on the full order linear evaluation model and 
demonstration in the nonlinear simulation. The next section 
documents how this flight test trajectory controller can be 
implemented for the zoom and pushover maneuver. 

5.1 TRAJECTORY MODELING 

The zoom and pushover trajectory has been designed to 
maximize the time at the target angle-of-attack, or alternatively 
maintain a constant vertical acceleration (with zero horizontal 
acceleration) . 

A parabolic flight path with apex at z Q has the property 

that 


1/2 

tan y = tan(0 - a) = £(z - z Q ) , (5.1) 

where y is the flight path angle and £ is a constant of pro- 
portionality, which defines the size of the parabola. An approxi- 
mate model can be written as 

1/2 

9 - a = £(z - Zq ) , ( 9-a) << 90° . (5.2) 

In addition, the following conditions must hold 

<f> = 0 . (5.3) 
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M and a are given at the apex. Our controller design will try 
to achieve the speed time history of a ballistic trajectory, such 
that the thrust nominally balances drag. With such a time history, 
the thrust would be automatically stabilized for the target angle- 
of-attack with desired Mach number and altitude at the apex of 
the parabola. Thus, the Mach number at any altitude is given by: 

M 2 “ Mq = 2g(z - z o )/c 2 , (5.4) 

where c is the speed of sound which is assumed constant on the 
trajectory (a variable value can also be accommodated). The 
equation can be linearized as 

2(M - M o )M q = 2| (z - z Q ) , (5.5) 


or 


M - M o ■ r^2 (z - V • 

M c 
o 

Note that the correct value of the speed will be reached at the 

apex because when z = z , M will be forced to M . 

o’ o 

£ is selected such that the lift defined by Mach number and 
angle-of-attack at the apex exactly balances gravity and the cen- 
trifugal force. A large value of £ will be needed when the 
target angle-of-attack or Mach number is small and vice versa 
(because small lift will require a tighter curve). It can be 
shown that with vertical acceleration a z and speed V at apex, 
£ can be selected as: 

V 2 £ 2 

■ L “2 i - = g - a z (5.6) 

2 

since the radius of curvature at the apex is 2 /£ . 
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Zoom and pushover is defined with 0-a and M as func- 
tions of altitude and the bank angle held at zero. 0-a and M 
time histories are two of several outputs of interest 

y = [a z , a, q, y, M] , 

which can be generated from a linear model as a function of time 
(see section 3). 

Since the thrust is stabilized, the exact final altitude 
cannot be achieved in the presence of gusts, hence we program 

y _ ( t ) and M (t) as functions of time rather than altitude. 

c C 1/2 
On a parabolic trajectory (z - z Q ) is linear in time. Since 

the trajectory command is open loop (due to the stabilized thrust 

constraint), the nonlinear relation between y(t) and M(t) can 

be used as shown below. The commanded output is then 


y c (t) = 


a Z c ( t) 

a c (t) 

q c (t) 

YcU) 
M c ( t > 


a c (a 2 c >Y c ,M,q ro ) 

0 


- Vtf + 

gR 0 2 
-2 2 [Y c (t)] 2 + M 0 


1/2 


V /p 6 - Y i 

(P 7 Y c ( t > + p 5 )1/2 


The angle-of-attack is not arbitrary since it is determined by 
the commanded a z and y and the actual Mach and dynamic 

C L/ 

pressure. The flight path angle is programmed from an initial 
value to zero at the apex. 

A summary of the zoom and pushover maneuver algorithm 
development is: 

1) Choose a , a desired apex angle-at-attack , and 
z Q , the desired apex altitude. 
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(2) Trim the aircraft at the apex altitude and angle 
of attack. 

(3) Find a final Mach number M q (which must be slower 
than the trimmed value to be on a parabola at the 
apex). This value cannot be arbitrary since the 
lift and drag are also functions of velocity. The 
force balance equations are: 



L - T sin a + W , 


0 = - D + T cos a 


We seek the apex velocity and thrust from 


+ c T ps m 
L , m 

o — + vr * + sin a 

2 R o 


- c D ps 


cos a 


W 


0 


2 

The acceleration V q /R o specifies the shape of the 
parabola and hence the constants in y (t) and 
M c (t). 

(4) Chose either the initial altitude or the time of 
flight (within an altitude range of approximately 
constant speed of sound, one specifies the other). 
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(5) 


( 6 ) 


Find the initial flight path angle and Mach number 
for the desired initial altitude, z^, from 


M. = 

l 


% <z - zp + 
c 


1 

tan y i = C(z Q - z i ) 2 

Regulate (a - a 2 ), (a - a (t)), q, (y - y (t)), 

o ° c 

and (M - M c (t)) about zero with a robust linear 

controller. 



5.2 ZOOM AND PUSHOVER (ZAPO) LINEAR MODELS 


This subsection describes linear models for: 1) the F-15 
aerodynamics, MCS/CAS, and sensor dynamics; and 2) the piloted 
aircraft with reduced order models for both. 

The trimmed condition of the F-15 used in the simulation 
SIMII is shown in Table 5-1. 


Altitude 

h 

o 

32,000 ft. 

Angle of Attack 

a o 

11° 

Pitch Angle 

e 

o 

11° 

Mach Number 

M 

o 

.4775 

Velocity 

V o 

470.9 ft/sec 

Thrust 

T o 

6270 lb f 

Weight 

W o 

40,700 lb f 

Stabilator 

«sb o 

-4.76° 

Dynamic Force 

qs 

55,821.1 lb f 


Table 5-1. ZAPO Trim Condition* 


♦Pitch, roll, and yaw CAS on, flaps up, with autopilots off. 




The six nondimen$ional aerodynamic coefficients and their 
associated dimension stability derivatives reflected into the 
appropriate body fixed axis are given below in Table 5-2. 


The MCS/CAS equations (ignoring servos at 30, 60, and 100 
radians/sec) from Figure 3-3 are 


<5gg - -,335gg + 15. 5q + 15. 5q 


^SB 2 5> 5 SB 2 + 1;L * 49 S e 


*SB, 


= .882q + . 4 a 


= .882q + 5.85 (a - q) 


SB 


■20. Esb + 20/57. 3(« SBi + <5^ + « SB 


+ . 6832q - 4.476 5 ) . 


The sensor and thrust equations are 
6 M = -.1<5M 


•M 

. i(«s\+ x (fZo\ 

\ M / ,x \ M / \ M / 
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ribriDIriENSIONAL CLjCHtCNjCIijClF i \ - 

-3.1M07E-13 - i ♦ 28433E-02 1..4240E-0B 1,10237 £-0j 7.05545E-01 4.44849E-12 


/.LfcO!'! CHBGDY CD CLFT CYBGDY 

Angle of attack 

36t07E-02 -0,l5745E+0i -v . 22784E-03 -0.61949E+02 -0.16598E+03 0.15015E-02 

Pitch rate 

G 0 , GOOvOE+OO — 0.91234E+00 0.00000E+00 O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+Ov 

Vertical velocity 
Vertical acceleration 

U -0.76929E-Q5 -0.32S19E-02 -0.47492E-06 -G.19277E+00 -0*3i496E+00 G.31297E-05 

hNZ 0* OOOOOE+OO -0.22116E-02 O.OOOOOE+OO O.OOOOOE+OO -0.30785E-01 O.OOOOOEtOO 

Angle of attack rate 

ALFIi O.OOOOOE+OO -0.76451E-01 O.OOOOOEtOO 0.00000E+00 O.OOOOOE+OO O.OOOOOE+OO 

Total velocity 

VELD -0 . 35605E-06 0.23469E-03 -0.35343E-07 -0.18301E-02 -0.40951E-02 0.50076E-06 

Velocity along A/C nose 

IJ 0.14954E-05 0.63794E-03 0.92315E-07 G.3/470E-01 0.61226E-01 -0.30835E-06 

Elevator 

DH o.OOOOOE+OO -0.33012E+01 -0.42944E-05 -0.12220E+02 -0.25475E+02 O.OOOOOE+OO 

Sideslip 

DEm -0.1O314E+O2 -0.30038E-06 0.1066QE+01 0.00000E+00 O.OOOOOE+OO -0.42589E+02 

Roll rate 

f -0.9143AE+00 0.00000E+00 -0.90524E-02 O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO 

Yaw rate 

R 0. 1-.7S1E+01 0.00000E+00 -0.25928E+00 O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO 

Lateral velocity 

V -0.83963E+01 -0.26201E-06 0.2ii48E+00 0.00G00E+00 O.OOOOOE+OO -0.37148E+02 

Aileron 

DO 0.25576E+01 0.74124E-09 0.10425E-01 0.00000E+00 O.OOOOOE+OO 0.10509E+00 

Rudder 

[ir — . 23370E+00 -0.58042E-Q7 -0.10299E+Q1 0.00000E+00 O.OOOOOE+OO -0.82293E+01 

Differential rudder 

J;T 0.33509E+01 -G.27126E-07 0.45272E+00 O.OOOOOE+OO 0.00000E+00 -0.38474E+01 


Table 5-2 


F-15 Aerodynamic Model 




The matrices of the linear model 


x = Fx + Gu + Tw 


y = H x , z = H X. + v 

y z 


are given below where the states 


x 


6V, 6a, q, 66, 6 gB , 





are velocity, angle-of-attack, pitch rate, pitch, four MCS/CAS 
states, lagged Mach number, and lagged thrust. The inputs 


u = 


6 e 6 T c/Mj ’ w = a gust ’ 

/ X" 4 2 

and w = N , 10, .1125 rad -sec 


are commanded elevator, thrust acceleration, and angle-of-attack 
gust. The outputs 


y = 



<5y» 



are acceleration, angle of attack, pitch rate, flight path angle 
(60 - 6a), and Mach number. The measurements 


z 



q, 



> 


are acceleration, angle-of-attack, pitch rate, pitch, and Mach 
number with 


V = 


N^O, di 


1.755 


lag 
\ -6 


O 3 \-6 o 

.2588 ft /sec , 7.62 rad -sec 


rad /sec, 


\ -4 2 \ 
1493 rad -sec, 2.5 


-6 


-sec 


The equations are scaled so that all angles, angular rates, and 
angular accelerations are in units of (.01) radians; velocities and 



Mach 


accelerations are in feet, with seconds the time unit, 
number is scaled in units of (.001) Mach. 


The linear system matrices are 


COLUMNS 

1 THRU 

8 


0.0019 

-0.3072 

0.0000 

-0.3220 

0.0009 

-0.3420 

1.0000 

0.0000 

0.0235 

-1.5480 

-0.9890 

0.0000 

0.0000 

0.0000 

1.0000 

0.0000 

0.3643 - 

•23,9940 

0.1705 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0156 

-3.3660 

-0.8723 

0.0000 

0.0056 

-0.3692 

-0.2359 

0.0000 

0.1014 

0.0000 

0.0000 

0.0000 

0,0000 

0.0000 

0.0000 

0.0000 

COLUMNS 

9 THRU 

10 


0.0000 

0.9816 



0.0000 

0.0019 



0.0000 

-0.0001 



0.0000 

0.0000 



0.0000 

0.0000 



0.0000 

0.0000 



0.0000 

0 , 000 ( 



0.0000 

0.0000 



0.1000 

0.0000 



0.0000 

-0.1000 




0.0000 

0.0000 

0.0000 

- 0.1222 

0.0000 

0 . 0000 

0,0000 

- 0.0541 

0.0000 

0 . 0000 

0.0000 

- 3.8000 

0.0000 

0.0000 

0 . 0 ' j 'j \j 

0 . 0000 

- 0.3330 

0.0000 

0.0000 

- 58,9000 

0.0000 

- 5.0000 

0.0000 

0,0000 

0,0000 

0.0000 

0.0000 

- 3,6681 

0.3490 

0.3490 

0.3490 

- 20.9063 

0.0000 

0.0000 

0,0000 

0.0000 

0.0000 

0.0000 

0,0000 

0.0000 


0,0000 0.0000 

0.0000 0.0000 

0.0000 0.0000 

0.0000 0.0000 

0.0000 0.0000 

11.4900 0.0000 

0.0000 0.0000 

-1.5620 0.0000 

0.0000 0.0000 

0.0000 0.1000 
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GM 


- 0.3072 
- 0.3420 
- 1.5480 
0.0000 
- 23.9940 
0.0000 
- 3.3660 
- 0,3692 
0.00 00 
0 . 0000 


COLUMNS 

1 THRU 

s 



- 0.0041 

- 1,6105 


0,0000 

0,0000 

0.0000 

1.0000 


0,0000 

0.0000 

0.0000 

0.0000 


1,0000 

0.0000 

0.0000 

- 1.0000 


0.0000 

1 .0000 

0.0000 

0.0000 


0.0000 

0,0000 

COLUMNS 

9 THRU 

10 



0.0000 

0.0000 




0.0000 

0.0000 




0.0000 

0.0000 




0 . 0000 

0.0000 




1.0000 

0.0000 





COLUMNS 

1 THRU 

8 


- 0.0041 

- 1.6105 

0 . 0000 

0.0000 

0 , 0000 

1.0000 

0.0000 

0.0000 

0,0000 

0.0000 

1 . 0000 

0.0000 

0.0000 

0.0000 

0.0000 

1,0000 

0.0000 

0.0000 

0,0000 

0.0000 

COLUMNS 

9 THRU 

10 


0 . 0000 

0.0000 



0.0000 

0.0000 



0 . 0000 

0.0000 



0.0000 

0,0000 



1.0000 

0.0000 




0.0000 

0.0000 

0,0000 

- 0,2548 

0,0000 

0.0000 

0 .0000 

0,0000 

0,0000 

0.0000 

0,0000 

0,0000 

0.0000 

0.0000 

0.0000 

0,0000 

0.0000 

0,0000 

0.0000 

0.0000 


0.0000 

0.0000 

0 . 0000 

- 0,2548 

0.0000 

0.0000 

0.0000 

0 . 0000 

0,0000 

0.0000 

0.0000 

0,0000 

0.0000 

0,0000 

0,0000 

0.0000 

0.0000 

0,0000 

0.0000 

0 , 0000 


4C 



Table 5-3 shows the eigenvalues of the aerodynamics alone 
as well as the aerodynamics as modified by the MCS/CAS. For 
level flight (the final condition at the apex of the ZAPO 
maneuver) at such a high angle-of-attack , the MCS/CAS (without 
the autopilot) makes the phugoid mode unstable. The time 
constant is on the order of the maneuver flight time, however, 
and the flight test trajectory controller (FTTC) will compensate 
it as well. 


Aerodynamics, sensor lag 
and MCS/CAS 

Aerodynamics alone 

0.0000 1 O.OOOOi 

-0.0042 + 0 . 0399 i 

0.0179 + 0.0259i 

-0.0042 - 0.0399i 

0.0178 •• 0.0259i 

-0.6622 1 1.19881 

-0.1000 t O.OOOOi 

-0.6622 - ] . 1988: 

-0.1000 + O.OOOOi 


-0.329S + 0.00001 


-1.2223 + 1.2490i 


-1.2223 - 1 . 24901 


-5.0000 + O.OOOOi 


-19.8335 + O.OOOOi 



Table 5-3. Eigenvalues for ZAPO Linear Model 


5.2.1 Nine State Aircraft Model 

Since the thrust is stabilized in the zoom and pushover 
model, the engine lag state, and thrust control will be dropped 
in the subsequent analyses. 

Frequency responses of the five outputs to sinusoidal 
elevator inputs are shown in Figures 5-1 and 5-2. 
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OUTPUT GAIN (see Sec. 
for respect, units) 


CM 

in 


1 . 

.1 

.01 

.001 

.0001 

.00001 

.000001 

.0000001 



.01 


Figure 5 - 1 , 


.1 1 . 

FREQUENCY (rad /sec) 


10 . 


Output Magnitude Responses to Elevator 
Aerodynamic and MCS/CAS Model 


Nino State 



MCS/CAS 
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IMAG PART (rad/sec) 


The a(s)/6 e (s) transfer function has the dip at low 
frequency because of zeroes produced by the MCS/CAS which 
attempt to mitigate the phugoid mode (see Figures 5-3 and 5-4) 


The high phase response is caused by the unstable phugoid 


mode . 


0.03 


0.02 


0.01 


0.00 


-.01 



REAL PART (rad/sec) 


Figure 5-3. a(s)/<$ e (s) — Low Frequency Poles and Zeros 

State Aircraft and MCS/CAS 



IMAG PART (rad/sec) 



REAL PART (rad/sec) 

Figure 5-4. a(s)/6 e (s) — Mid-Frequency Poles/Zeros — Nine State 

Aircraft and MCS/CAS 
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The residues of these transfer functions are given below in 
Table 5-4. 


Eigenvalues 

a 76 
Z e 

a/5 

e 

q/5 e 

Y/« e 

M/6 

V 

0.0000 

+ 

O.OOOOi 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0178 

+ 

0 , 0259i 

-0.0026 

0.0018 

0.0006 

0.0150 

0.0183 

0.0178 

- 

0.02591 

0.0009 

0.0002 

-0.0002 

0.0110 

-0.1753 

-0.1000 

+ 

O.OOOOi 

0.0000 

0.0000 

0,0000 

0.0000 

-0.0163 

-0.3298 

+ 

O.OOOOi 

-0.0003 

0.0002 

0.0000 

-0.0002 

-0.0001 

-1.2223 

+ 

1 . 2490i 

-0.0800 

0.0441 

0.0879 

-0.0203 

-0.0020 

-1.2223 

- 

1 .2490i 

-0.1557 

0.0999 

-0,1441 

-0.0063 

0.0010 

-5.0000 

+ 

O.OOOOi 

0.0104 

-0.0607 

0.3012 

0.0004 

0.0000 

-19.8335 

+ 

0.00001 

0.4704 

0,0146 

-0.3897 

0.0050 

0.0001 


Table 5-4 Residues of Elevator Transfer Functions 
(Nine State Aerodynamics and MCS/CAS) 


These residues indicate that a good model can be achieved 
with five states, by including the phugoid and short period modes 
and the control system lag at 5 radians/sec. Such a model would 
be useful for flight test trajectory control signals fed directly 
into the aircraft control system through the autopilot. This 
reduced order model is detailed in the next subsection. 


5.2.2 Five State Aircraft Reduced-Order Model 

The nine state model of the previous subsection, reduced 
to match at zero frequency by the method of Section 4.1 gives 
the following equations 


x v. = F x + Gu + w , 
r r r r r 


y = H y x r + J u u » 
* r r 


z — H x + J w, 
z r w ’ 
r r 


with the matrices given by 
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0.0178 

0.0259 

0.0000 

0.0000 

0,0000 

-0.0259 

0.0178 

0.0000 

0.0000 

0.0000 

0,0000 

0.0000 

-1,2223 

1.2490 

0.0000 

0.0000 

0,0000 

-1,2490 

-1.2223 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-5.0000 


0.2905 
-0,0541 
G r -0.1374 

0.0904 
12.1783 

0.5576 

0.1872 

2.0404 

-0.7345 

0.0000 


-0.0092 

-0.0014 

-0.1140 

-1.0581 

0.0009 

0.0060 

-0.0017 

0.1101 

0.6546 

-0.0050 

0.0021 

0.0004 

-0.9280 

-0.4379 

0.0247 

0.0431 

-0.0459 

0.0822 

-0.0998 

0.0000 

0.1693 

0.5720 

0.0133 

-0,0014 

0.0000 



-0.0696 

0.0012 

-0.0197 

- 0.0002 

-0.1611 


-0.0092 

-0,0014 

-0.1140 

-1.0581 

0.0009 

0.0060 

-0,0017 

0.1101 

0.6546 

-0.0050 

0.0021 

0.0004 

-0.9280 

-0,4379 

0.0247 

0.0371 

-0,0442 

-0.0279 

-0.7544 

0.0050 

0.1693 

0.5720 

0,0133 

-0.0014 

0.0000 


0.4187 

-1.8776 

0.0034 

2,1376 

-2,5170 


The frequency response of this five state model is shown in 
Figure 5-5. The model is only good up to about 4-12 Radians/sec 
and only up to .1 radians in the mach response which is primarily 
due to its own lag at .1 second which was not included in the 
model . 
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eg 



Figure 5-5 Transfer Function Magnitudes - Five State Aircraft 
and MCS/CAS (Reduced Order Model Matched at Zero 
Frequency) 
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5.2.3 Thirteen State Piloted Aircraft Model 


The four state pilot model of Section 3.5 is for a 1/s 2 
plant. The pilot adapts his own gain to make the plant look like 
the dynamics he wants. At this point we can ignore the gain 
question of the pilot until after the compensator has been 
designed, then adjust the compensator gain to simulate the 
pilot's action. 


The unreduced design model is given by 


- 


— 

— 


“ — 


- 


- 

— 


- 

X A/C 

• 

= 

f a/c 

g a/c h p 


X A/C 

+ 

0 

+ 

r 

0 


“gust 

x„ 


0 

F„ 


x„ 


G 


0 

r 


w 

P 



P J 


L p J 


L p J 



p 


p J 


y = [H 

0] , 

y A/C 


z = [H 

0] + v . 

Z A/C 



The magnitude and phase responses of the five outputs to dis- 
played elevator inputs are given by Figure 5-6 and 5-7. These 
frequency responses are not surprising when looking at transfer 
function of the linear pilot model (see Figures 5-8 and 5-9). 

The response of this thirteenth order system to gust inputs 
is shown in Figure 5-10. 

The residues of the elevator transfer functions are the 
basis for model reduction. They are shown in Table 5-5. 

The pilot model dominates, so these 4 modes together with 
the unstable phugoid mode comprise the reduced order design 
model described in the next subsection. 
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eg 

in 



FREQUENCY (rad /sec) 

Figure 5-6 Response to Displayed Elevator Command - 13 State 
A/C and Pilot Model 



Figure 5-7 Phase Response to Displayed Elevator Command - 
13 State A/C and Pilot Model 
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)T PHASE (de 






OUTPUT GAIN (see 


CM 



Table 5-5 Elevator Residues - Thirteen State Piloted Aircraft 


Eigenvalues 


a /6 
z e 


cx/6 

e 


q/5 e 


Y/6 e M/(S e 


0.0000 

+ 

O.OOOOi 

-3.7851D-32 

-3, 1477D-17 

-8 . 6754D-31 

0.0178 

+ 

0.0259i 

6 , 7004D-06 

-4 , 8157D-06 

-1 , 5159D-06 

0.0178 

- 

0.0259i 

-2 . 6021D-06 

-2 . 7827D-07 

6 , 7240D-07 

-0.1000 

+ 

O.OOOOi 

-3 .82060-21 

2 . 5099D-21 

2 , 0833D-21 

-0.3298 

+ 

O.OOOOi 

3. 2125H-07 

-2 . 1251D-07 

1.8614D-09 

-1.2223 

+ 

1 ,2490i 

1 .24718-03 

-8 . 0742D-04 

1.3152D-03 

- 1 .2223 

- 

1.2490i 

-8.3688D-04 

4 . 7492D-04 

6 ♦ 0570D-04 

-5.0000 

+ 

O.OOOOi 

-2 , 7288D-04 

1 , 5872D-03 

-7 . 8779D-03 

-8.0000 

+ 

0.00001 

-1 , 0129D+00 

-4.0442D-01 

3. 4505H+00 

-9.5240 

+ 

O.OOOOi 

1 , 2567D+01 

2 . 7990D+00 

-2 . 9326D+01 

-10.0000 

+ 

0.00001 

-1 , 3282D+01 

-2 . 5525D+00 

2.8345D+01 

-14.0000 

+ 

O.OOOOi 

2 . 2990D+00 

1 . 7493D-01 

-2 ♦ 9372D+00 

-19.8335 

+ 

O.OOOOi 

-5 . 7243D-01 

-1.7783D-02 

4 , 7425D-01 


-1 . 179BD- 17 
-4 ,01 16D-05 
-2.7340D-05 
2 • 7531D-20 
2 . 0687D-07 
3.331 ID-05 
-1 * 7944D-Q4 
-1 . 1590D-05 
-2.6886D-02 
2 . 8 0 2 1 D - 0 1 
-2.8205D-01 
3.4872D-02 
-6,12910-03 


8,47300-15 
-3.19570-05 
4.58710-04 
3.55780-05 
7.62950-08 
-1,05490-05 
-1.56580-05 
4 , 15040-07 
-8.58580-04 
7.16130-03 
-6 . 80690-03 
5 .81820-04 
-7. 12870-05 
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5.2.4 Six State Piloted Aircraft Reduced-Order Model 


The reduced piloted aircraft model is of the same form as 
the equations of Section 5.2.2, with the following matrices 


o 

o 

^>1 

00 

0,0259 

0,0000 

0.0000 

0.0000 

0.0000 

-0.0259 

0.0178 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-8.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-9,5240 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0,0000 

-10.0000 

0,0000 

0,0000 

0.0000 

0,0000 

0.0000 

0.0000 

-14.0000 


1.0D+04 * 


G = 

0,0000 

0,0000 

0.0985 

-4.2976 

-0.1120 

0.3172 






r = 

0.5576 

0.1872 

0.0000 

0.0000 

0.0000 

0.0000 

-16.2773 

2.9859 

0.0000 

0.0000 

44.2543 

0.0000 





H y ■ 

"0.0092 

0.0060 

0.0021 

0.0431 

0.1693 

-0.0014 

-0.0017 

0.0004 

-0.0459 

0.5720 

0.0577 

0.0230 

-0.1965 

0,0015 

0.0000 

0.0164 

0,0037 

-0.0383 

0,0004 

0.0000 

-0.6655 

-0,1279 

1.4202 

-0.0141 

-0.0003 

-0.0407 

-0.0031 

0.0520 

-0.0006 

0.0000 

J = 
u 

1.6135 

0.0398 

-1.2966 

0,0208 

-0.0192 






H z = 

-0.0092 

0.0060 

0,0021 

0.0491 

0.1693 

-0,0014 

-0,0017 

0,0004 

-0.0476 

0.5720 

0,0577 

0,0230 

-0,1965 

0,0246 

0,0000 

0,0164 

0,0037 

-0,0383 

0.0040 

0.0000 

-0.6655 

-0,1279 

1.4202 

-0,1420 

-0.0003 

-0.0407 

-0,0031 

0.0520 

-0,0037 

0,0000 

Ch 

S: 

II 

1.5377 

-0.9490 

0.0185 

-0,5339 

-2,5086 

6.8958 

-2.0947 

-7.2272 

-1.3389 

9,2763 
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OUTPUT GAIN (see Sec. 
for respect, units) 


The frequency magnitude response to elevator inputs for this 
sixth order reduced order model is given in Figure 5—11. This 
model is considerably different from the linear evaluation 
(thirteenth order) model above one radian per second. However, 
the model information at low frequency is quite good and, with 
frequency shaped control methods, can be used to develop a robust 
flight test trajectory controller. Reducing the thirteenth order 
model at a higher frequency does not significantly improve the 
flight path angle response which can be observed in Figures 5-12 
and 5-13. Figures 5-14 and 5-15 show the response of this reduced 
order design model to gusts and pilot motor noise. 



Figure 5-11 Magnitude Response to Displayed Elevator Command - 
Six State Reduced Piloted A/C Model (Reduced at 
Zero Frequency) 
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Lndino 


FREQUENCY (rad /sec) 

Figure 5-12 Response to Displayed Elevator - Six State Piloted 
A/C Model (Reduced at 1 Rad/Sec) 



FREQUENCY (rad/sec) 

Figure 5-13 Response to Displayed Elevator - Six State Piloted 
A/C Model (Reduced at 10 Rad/Sec) 




OUTPUT GAIN fin respective units) OUTPUT GAIN (in respective 
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Figure 5-14 Response to Gust - Six State Piloted A/C Model 
(Reduced at Zero Frequency) 



Figure 5-15 Response to Pilot Motor Noise - Six State Piloted 
A/C Model (Reduced at Zero Frequency) 
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5.3 FLIGHT TEST TRAJECTORY CONTROL ( FTTC ) COMPENSATOR 

Using the reduced order design model of the piloted aircraft 
from the previous section, a standard linear quadratic regulator 
could be designed with the following quadratic costs on outputs 
and inputs, 

A y = diag [.01, .01, .01, 10, 1] ,B 6e = 1 . 

This yields the following control gains and closed loop regulator 
eigenvalues shown below, 

C = [1.4068 -1.7544 .0011 .0003 -.0106 -.0006] , 

with X ± (F r -G r C) 

-12.511 ± 1. 8195i 
-8.0244 ± . 9075i 
-.0210 ± . 0280i . 

The closed loop regulator response to incremental elevator is 
shown in Figure 5-16. 

The poor knowledge of the plant at high frequency as shown 
in Figure 5-16 cannot be overcome with regulator design. The 
next subsection shows a frequency shaped controller which uses 
the high fidelity of the model at low frequency, gives acceptable 
performance at mid-frequency where the pilot works, and rolls off 
the controller at high frequency. 

5.3.1 Regulator Design 

A more robust controller with respect to the reduced order 
modeling errors can be obtained by using frequency shaped control 
concepts as described in Appendix A. The following performance 
index 


64 



OUTPUT GAIN (in respective units) 



FREQUENCY (rad /sec) 


Figure 5-16. Response to Elevator with Full State Feedback — Six 

State Piloted A/C 
(Plant reduced at zero frequency) 
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w 

lim 
w->°° -w 


W 

/ 


‘ + (,01 + w 2 ) B {e 


(1+W~) 


dw , 


can be implemented by the following equations 



X 


F 

O 

0 


O 

d_ 

dt 

Z 1 

= 

GH 

a 

z 

F 1 

G. J 
1 u 

a 

z 

dx + 

O 


u 


O 

O 

1 

<M 

Pm 


■ 

o 

E\3 

I 


u 


where F^ = 

-1 1 

’ G 1 = 

0 


-1 


1 


Fg • 1 > Gg .1 , 


and F, G, H a , and Ju & are given in Subsection 5.2.4, 
z z 

with the performance index as 

T , 


lim gTj, 

o 


~ r 

2T J 


2 ^ - 2 \ 
z 1 + u I dt 

^ 1 1 


The control law is given by 


u = -C 1 x - CgZ-^ - CgU , 


with feedback gains and closed loop regulator eigenvalues as 
follows : 

COLUMWS 1 THRU 3 

C A = f C l ,C 2 ,C 3-' = 14 ♦ 2516 -14.5265 0,0000 0.0000 0.0000 0,0000 0,0009 

COLUMNS 9 THRU 9 
0,7124 


0.0003 
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OUTPUT GAIN (in respective units) 


X j .(F A -G A C A ) = -14.0000 - G.OOOOi 

1 ' " -9.5240 - 0,00001 

-10.0000 - O.OOOOi 
-8.0000 - O.OOOOi 
-1.0000 + 0.0032i 
-1.0000 - 0.00321 
-0.0671 + 0.252*1 
-0.0671 - 0. 252<?i 
-0.0014 + 0.0000: 


The control law for u can be obtained from the ninth state 
of the compensator as shown in Subsection 5.3.3. Our control 
u = 6 e is primarily used to control flight path angle. By 
rolling off the acceleration, the derivative of flight path 
angle, at high frequency, the controller achieves good roll off 
in all of the outputs as shown in Figure 5-17. 



Figure 5-17. Closed Loop Regulator Response to Gust 
(Six piloted A/C states with three shaping states) 
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5.3.2 Estimator Design 


An estimator can be designed for the six state piloted air- 
craft model with the gust and measurement power spectral densities 
given in Subsection 5.2.3. Using all five measurements and 
w(t) = a st input the performance index is 

J = E [w T (t )Q - ‘*'w(t ) + v T ( t )R~lv( t )J , 

yielding an estimator feedback gain matrix K to implement the 
filter equations 


x = Fx + Gu + K(z-H 


z 


X 


J u) • 
u 


The gains and filter closed loop eigenvalues are shown below. 


K 


0.1409 

-0.2728 

0.2119 

0,5487 

-0.0320 

0.0486 

-0.1244 

-0.0345 

-0,3664 

0,1011 

0.0000 

0 . 0000 

0 . 0000 

0,0000 

0,0000 

0.0000 

0.0000 

0.0000 

0,0000 

0,0000 

0.0000 

0,0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0,0000 

0,0000 


A.(F-KH ) = 
l z 

-0.0294 t 0.0295i 
-0.0294 - 0.0295i 
-14.0000 + O.OOOOi 
-8.0000 + O.OOOOi 
-9.5240 + 0.00001 
-10.0000 + O.OOOOi 
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5.3.3 Reduced Order Compensator 


The full state compensator is given by 



X 


F-kH 

z 

0 

G-kJu H 
z c 


X 


k 

d 

dt 

Z 1 

= 

GH a 

z 

F 1 

G l Ju a 

z 


Z 1 

+ 

0 


u 


-Ct 2 C i 

- G 2 C 2 

F 2 _G 2 C 3 


u 


0 


with H =[001]. 
c 

The feedforward of the control into the measurements requires 
a negative feedback in the compensator. The residues of this 
transfer function are given below in Table 5-6 


Eigenvalues 

U £/a 

z 

°D/a 

U D/q 

U D/ Y 

U D/M 

-0.0793 t O.OOOOi 

-0.2075 

0.5248 

0.1258 

1.4533 

-0.4092 

-0.0807 9 0.22891 

0.2074 

-0.5247 

-0.1259 

-1,4536 

0.4093 

-0.0807 - 0.22891 

-0.5675 

0.9060 

-1.5388 

-5.7511 

0.8434 

-1.0000 t 0.0033i 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-1.0000 - 0 . 0033 l 

0,0005 

-0,0007 

0.0014 

0.0053 

-0.0008 

-8.371? 9 0.65901 

-0,0005 

0.0008 

-0,0013 

-0.0049 

0.0007 

-3.3710 - 0.6590i 

0,0019 

-0,0030 

0.0052 

0.0193 

-0.0028 

TO. 9263 9 O.OOOOi 

0.0006 

-0.0010 

0,0017 

0.0065 

-0.0010 

-12.8433 9 O.OOOOi 

-0,0001 

0.0002 

-0,0003 

-0.0012 

0.0002 


Table 5-6. Full State Compensator Residues. 


*U D is displayed control (elevator) command 
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This compensator can be well approximated with a three state 
compensator using modal reduction just as the original model was 
reduced (see Section 4.1). The reduced order compensator is 
given by 


x 

c 


F x 
c c 


+ G z 
c 


u = H x , 
c c * 


with 


FC 


-1.9338B-01 

2.5536B-01 

-2.3696B-01 

2.0763B-01 

-5.5726B-01 

3 . 0108D-01 

GC 


-8.9338B-01 

1.5645B+G0 

-1 » 1924D+00 

2.2681B+0G 

-5.9027D-01 

3.0340B-01 

HC - 



1.3422D-01 
6.4751B-G2 
-2. 5491H-01 


-1.9311B+00 -6.5153B+00 8.1542D-01 
-1.9376D+G0 -5.3908B+00 4.2135B-01 
-2.0951D+00 -8.5362B+Q0 1.39270+00 


4. 4421D-01 -2.4000B-01 -1.8743D-01 , 


The magnitude and frequency responses of these five transfer 
functions are shown in Figures 5-18 and 5-19. On the same scales 

the ninth order frequency response plots are indistinguishable as 
would be expected from the size of the residues. 


5.3.4 Controlled Plant Dynamics 

The compensator can be further analyzed by investigating 
the effect of the dynamic compensator on the full thirteen state 
aircraft and MCS/CAS model, G(s). Viewing the pilot as a 
compensator, P(s), for this system, we are interested in the 
transfer function from the aircraft stick to the displayed 
elevator command to which he responds (see Figure 5-20 below). 
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OUTPUT PHASE (deg) OUTPUT GAIN (in respective units) 



.01 .1 1 . 10 . 

FREQUENCY (rad /sec) 

Figure 5-18. Throe State Reduced Order Compensator 

Magnitude Responses. 


FREQUENCY (rad /sec) 

Figure 5-19. Three State Reduced Order Compensator 

Phase Responses. 
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Aircraft 


FTTC 



Figure 5-20. Controlled Plant Dynamics 


The state equations for the controlled plant are 



The frequency response of the controlled plant is shown in 

Figures 5-21 and 5-22. In the region of .5-2 radians/sec the 
2 

plant is k n /s and in the region of 2.-20 radians/sec it is 
1 3 3 

approximately k„/s . A k/s plant is clearly harder to 

^ 2 

control than a k/s plant, however the plant gain at the 
upper range of these frequencies is so low (shown on the log-log 
scale here) that the change in dynamics is not significant. The 
pilot model is based on negative feedback, hence the phase of the 

plant should ideally remain close to zero, and as shown in 
s 

Figure 5-22 remains between -20° to -50° in the .5-2 radian/sec 
frequency range. 
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. CONTROLLED PLANT GAIN 



Figure 5-21. Controlled Plant Magnitude Response 
(9 State A/C + 3 State Compensator) 



Figure 5-22. Controlled Plant Phase Response 
(9 State A/C + 3 State Compensator) 
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5.4 LINEAR CONTROL LAW EVALUATION 


Before evaluating the compensator in the nonlinear simulation 
the compensator must meet linear evaluation tests for stability, 
disturbance rejection, and other desirable controller properties. 

The effect of sampling rates can be evaluated much more inexpen- 
sively on the linear model than with extensive nonlinear simulation. 


5.4.1 Stability 

Stability evaluation of the compensator can be systematically 
checked by computing the eigenvalues of 

1) the full-state compensator and design model (these should 
be the regulator and estimator eigenvalues) 

2) the full-state compensator and evaluation model 

3) the reduced order compensator and the design model, and 

4) the reduced order compensator and the evaluation model. 

In each case an augmented dynamics matrix is formed 



G H 


c z 



F 


x 


J 


and the eigenvalues computed. The eigenvalues for the second and 
final tests are shown in Table 5-7 below. 


5.4.2 Disturbance Rejection 

The performance of the regulator alone assuming full state 
feedback in rejecting gust disturbances can be evaluated with the 
following equations 

(F-GC)X + X(F-GC) ' = - r Q gust rT 

Y = H XH . 

y y 
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22nd Order 
Evaluation 
Model 

19.8241 + O.OOOOi 
•14.0000 + O.OOOOi 
-5.0023 - O.OOOOi 
-1.2176 + 1.2470i 
-1.2176 - 1.2470i 
-1,0000 - 0.0033i 
-1.0000 + 0.0033i 
-0,0666 - 0.2403i 
-0.0666 + 0.2403i 
-0.3298 - O.OOOOi 
-0,1208 - 0.00001 
-0.0272 - 0.027U 
-0.0272 + 0.0271i 
-0.0024 + O.OOOOi 
•14.0000 - O.OOOOi 
10.0013 t O.OOOOi 
-9.5222 - O.OOOOi 
-8.0006 t O.OOOOi 
■10.0000 + O.OOOOi 
-9.5240 - O.OOOOi 
-8.0000 - O.OOOOi 
0.0000 + O.OOOOi 


16th Order 
Evaluation 
Model 

-19,8241 + O.OOOOi 

-14.1338 - O.OOOOi 
-9.8726 I 1.0553i 
-9.8726 - 1.0553i 
-7.6555 + O.OOOOi 
-5,0024 - O.OOOOi 
-1.2176 I 1.2470i 
-1.2176 - 1.2470i 
-0,0668 + 0,2410i 
-0.0668 - 0.2410i 
-0,3298 - O.OOOOi 
-0.1206 f O.OOOOi 
-0,0270 + 0.0271i 
-0,0270 - 0.0271i 
-0.0025 + O.OOOOi 
0.0000 + O.OOOOi 


Table 5-7 


Closed Loop Eigenvalues 



The rms values are the square roots of the diagonals of Y, 

2 

shown below with a gust spectral density of .0369 deg -sec: 

YRHS ; 

1.0737D-01 g 
1.3364D-01 deg 
4,83790+00 deg /sec 
8. 96510-02 deg 
1,09810+00 . 001 Mach 

URNS - 

6, 50740-02 inches (Longitudinal stick motion) 

This control law has not been optimized but can be improved 
with additional design iterations. 

5.4.3 Sampling Rate Effects 

The effect of sampling rates can be evaluated by discretizing 
the aircraft model as well as the pilot model and FTTC compensator, 
at the least common multiple of all the sampling rates involved. 

The dynamics augmented with two delays for the display update and 
fifteen delays for the sensor update and the appropriate discrete 
feedback compensation, can then be evaluated for average behavior 
by solving the discrete form of the symmetric Lyapunov equation. 

The discrete model for the compensator (1/A = 160 Hz) is, 
x c (k+l) = F D x c (k) + G D z(k) 
u(k) = H D x c (k) , 

where these matrices are given by 
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FD * 

9.98790-01 1.59690-03 8.3799D-04 

-1.48180-03 1. 0013H+00 4. 04020-04 

-3,47940-03 1,87370-03 9.98410-01 

GO - 

-5.58770-03 9.78590-03 -1.20770-02 -4,07450-02 5.09910-03 

-7,45380-03 1,41780-02 -1.21120-02 -3.36950-02 2.63320-03 

-3.68360-03 5.01360-03 -1.30740-02 -5,32700-02 8.69100-03 

HO = 

4.44210-01 -2.40000-01 -1.87430-01 . 

The eigenvalues in the z-plane of the compensator are: 

9.99490-01 +1.43430-031 
z. = 9.99490-01 -1.43430-031 
1 9,99510-01 -3.82010-211 . 


The discrete model for the pilot (1/A also 160 Hz) is 


9.39410-01 -2,16280+00 -6.82230-03 -1.36130-05 

0,00000+00 9,99960-01 6.23750-03 1.82900-05 

0.00000+00 -1.95100-02 9,94080-01 5.66090-03 

0.00000+00 -6. 03840+00 -1.83960+00 8. 15620-01 

GO = 

5.40520-03 

5.85150-04 

-1.25200-02 

2.23700-01 

HO - 

-5,61030+01 0.00000+00 0.00000+00 0.00000+00 

with z-plane eigenvalues, 

9.39410-01 +0.00000+001 
z = 9.51230-01 -3.31480-17! 

i 9,42210-01 +5.6503D-17i 
9.16220-01 -1. 94400-17! 
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5.5 NONLINEAR SIMULATION CONTROL LAW EVALUATION 


After linear evaluations, the control law can be more signi- 
ficantly evaluated by demonstration in the nonlinear simulation. 
This required closing the control loop of the simulation; the 
controls no longer need to be open-loop commands, but rather 
can be feedback from the simulation outputs, corrupted by a 
desired noise level. The details of how this FTTC simulator 
can be invoked are given in the next section. 

The zoom and pushover control law developed in this section 
was implemented on the modified version of NASA Dryden's SIMII 
simulation. Following the algorithm development steps outlined 
in subsection 5.1: 

Steps (1) — (2) : The F-15 trim values were chosen for the 

desired apex conditions 

a Q = 11 degrees 
z Q = 32,000 feet 

o 

Steps (3) — (4) : Assume an acceleration of V 0 / R 0 = 1° feet 

per second to climb a thousand feet in approximately 
15 seconds. Thus, the radius of curvature is approxi- 
mately 20,000 feet (for a that is 90% of the trimmed 
value ) . 

The relations in Section 5.1 give V Q = 406 feet /second and an 
initial flight path angle of = 17.55°. This is a little 

too far away from the linearization of the aerodynamics at 11°. 
Fixing y^ at 11 degrees gives the trajectory parameters shown 
in Table 5.8. The linear model trim condition is approximately 
half way in the maneuver velocity, and hence dynamic pressure 
range . 
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TABLE 5-8 ZOOM AND PUSHOVER DEMONSTRATION 
TRAJECTORY PARAMETERS 


z 

o 

Z i 

1000 feet 

R 

o 


52,933. feet 

£ 


6.147 10 - 3 

V o 


442. feet/sec 

M 

o 


.4491 

T 

o 


5544.6 lb f 

(g ~ 

V 

3.706 feet /second^ 

t o - 


23.23 seconds 

Y i 


11. ° 

V. 

l 


510.4 feet/second 

M . 

1 

.5175 

l 




To simplify the a_(t) command we keep it constant while 

c 

decreasing the flight path angle linearity in time under a stabi- 
lized engine power lever angle condition. The thrust itself is 
not stabilized because.it is also a function of Mach, as the 
aerodynamics are. This maneuver is an acceleration initial 
condition response, shown in Figure 5-28 through 5-32. 

The loss of velocity (see Figure 5-32) causes a higher angle- 
of-attack and acceleration seen in Figure 5-28 and 5-29. The 
velocity is lost because it is difficult to trim the aircraft in 
a condition where the zoom and pushover can begin. The ballooning 
effect, seen in the flight path angle response together with the 
variation of the thrust with velocity accounts for the loss of 
energy. It reaches a stable angle-of-attack condition. Initial 
trim would provide near perfect response. 

Use of throttle would give considerably better control, and 
would also enable modeling the trajectory commands as a linear 
system so that all the control loops could be closed. Random 
walk states could also be appended to model the effects of changing 
dynamic pressure and thrust. MIMO servomechanism control could 
also be used since the response is influenced by the trajectory 
commands significantly. 
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22 . 
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20 . 

15. 

18 . 

17 . 

16 . 

15 . 

14 . 

13 . 

0 . 3 . 6 . 9 . 12 . 15 . 18 . 21 . 24 . 

TIME(SEC) 

Figure 5-25 Pitch - Degrees (a Initial Condition) 

z 
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0 . 3 . 6 . 9 . 12 . 15 . 18 . 21 . 24 . 


TIME(SEC) 


Figure 5-26 Flight Path Angle - Degrees (a Initial Condition) 

z 


82 




Figure 5-27 Mach (a Initial Condition) 

z 
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5.6 ZOOM AND PUSHOVER ALGORITHM SUMMARY 


The development of a performance index and trajectory 
commands with the necessary initial for the zoom and pushover 
maneuver was described in trajectory modeling section. Reduced 
order models for the aircraft and piloted aircraft were described 
and a linear compensator synthesized to regulate about the 
necessary trajectory commands. This compensator was reduced 
to three modes, and its stability and disturbance rejection 
properties evaluated with the full order linear model (thirteen 
states) . 

The compensator and pilot models were discretized and imple- 
mented in the nonlinear simulation as described in the next 
section and the FTTC User's Guide (Appendix C). The initial 
condition response illustrates the desirability of closing the 
loop on the trajectory commands as well. This demonstration 
used an open loop trajectory command generator, since the thrust 
was stabilized. 

The particular control law, the linear compensator and the 
chosen trajectory command generators, is not optimized but 
illustrates all the aspects of the systematic design approach 
including the maneuver modeling, linear model development, linear 
control design, linear evaluation, and nonlinear evaluation with 
simulat ion . 
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SECTION 6 


ZOOM AND PUSHOVER ALGORITHM DOCUMENTATION 


The zoom and pushover maneuver FTTC is shown in block 
diagram form in Figure 6-1. The inputs and outputs used by each 
block are given along with references to sections of this report 
which document the equations implemented by each block. 

For the zoom and pushover with a sensor data rate of 10 Hz 
and display update rate of 53-1/3 Hz, the simulation can be inte- 
grated at the least common multiple of the sampling rates, in 
this case 160 Hz. The discrete compensator equations are updated 
every sixteen aircraft integration steps and the display is up- 
dated every three steps. (Note that the multiple of integration 
time step fcr the sensors and displays can be input to accommodate 
general multirate systems.) 

The compensator can have arbitrary structure including dyna- 
mics and feed forward gains. The pilot model is actually four 
independent models which affect the four pilot controls: elevator, 

aileron, rudder, and throttle. Any or all of these can be acti- 
vated and can have identical or different structure and/or param- 
eters which are defined with input data (see Appendix C). 

Neglecting for a moment the dynamics of the compensator, the 
corresponding regulator control law is of the general form 

u = - C x (x - x T - x c ) + u T + u c , 

where x m are the trim states and x are the commanded states 
from the trajectory generator. The overall flight test trajectory 
controller is a linear feedback compensator, with set point control. 
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00 

o> 


(Section 5.2) 



Figure 6-1 Zoom and Pushover Flight Test Trajectory Controller (FTTC) 
Block Diagram 









and an open loop command generator. If we generated the trajec- 
tory commands as outputs of a linear system, these loops could 
bo closed as well in a MIMO servomechanism. On the zoom and 
pushover maneuver, however, it is the constraint of stabilized 
thrust that requires the open loop trajectory commands. 

This FTTC structure can be used for general simulation tasks 
as well. For example, the pilot model could be bypassed with 
pure feed forward unity gain, the sensors disconnected, and 
pure feed forward on the compensator to use the trajectory gener- 
ator to input an open loop control function, such as a swept 
sine wave. 

The flowchart for demonstrating the FTTC in the zoom and 
pushover maneuver, and any other maneuver, in a modified version 
ol the SIMII simulation, is shown in Figure 6-2. The flowchart 
shows the sequence of operations in controlling a flight test 
maneuver and how the FTTC is interfaced with the NASA DRC 
simulation SIMII. The subroutine FTTC integrates the functions 
described in the block diagram (Figure 6-1). FTTC is called 
by the subroutine CNTRL when flight test trajectory simulation 
is requested, otherwise CNTRL implements an open loop control 
law, entered as segments of a time history. Only the subroutine 
associated with a particular sensor configuration (SENSOR) and 
a particular trajectory generator (TRAJGN) need to be recompiled 
for a different maneuver. All the others have a general structure 
suitable for all trajectories. These subroutine listings are 
given in Appendix D. 

The outputs from the various blocks of the FTTC are shown 
in Table 1 of Appendix C. These user-accessible variables can 
be output from the simulation with other time history variables 
of interest. 
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Aircraft equations 


Flight Test Trajectory Controller Flowchart for the Batch Simulation (SIMII) 

Demonstration 




SECTION 7 


SUMMARY 


Recent extensions to Linear-Quadratic-Gaussian (LQG) 
synthesis theory (Appendix A) applied to meaningful linear models 
(Section 5.2) with sufficiently flexible software tools (Appen- 
dix B) provide powerful techniques for designing Flight Test 
Trajectory Controllers (FTTC). The effort reported here has 
adapted and developed the necessary software and models to 
implement the zoom and pushover trajectory controller in the 
NASA-BRC batch simulation SIMII. The FTTC structure described 
in this report is sufficiently general to be used for other 
maneuvers, with the aid of' the FTTC User's Guide [Appendix Cl. 

The zoom and pushover maneuver controller design (Section 5) 
provides a detailed roadmap for the design procedure described 
in Section 4. This procedure involves a significant yet reason- 
able amount of engineering effort given sufficiently powerful 
software tools. Subsequent paragraphs summarize conclusions on 
the FTTC synthesis procedure, the linear models, the control 
design on the zoom and pushover maneuver, technical problems 
encountered and recommendations for further research. 

7.1 SYSTEMATIC METHODS FOR FTTC DESIGN 

The principal steps for systematic development of flight 
trajectory controllers can be summarized as planning, modeling, 
designing, and validating a trajectory controller. An important 
part of the planning is the trajectory specification. Major 
classes include single-vehicle and multiple-vehicle trajectories. 
Subclasses include constraints on combinations of load, speed 
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and altitude. General models suitable for trajectory control 
were demonstrated for the major subsystems, aerodynamics, kine- 
matics, engines, sensors and pilot, by linearizing about the 
flight path. The control-design procedure consists of model 
reduction of the linear models where possible, design of a con- 
trol law with a regulator /estimator or output feedback controller, 
and reduction of the controller complexity where possible for 
more robustness. 

7.2 LINEAR MODEL DEVELOPMENT 

Having access to a complete set of dimensional stability 
derivatives based on a detailed nonlinear F-15 aerodynamic model 
at arbitrary flight conditions is highly desirable. The specifi- 
cation of flight condition cannot be made completely arbitrary, 
however, because of the requirement to trim the aircraft, dis- 
cussed further below. These dimensional stability derivatives 
are projected onto both stability and body axes. The longitudinal 
MCS/CAS, reduced to a fourth order model, is dependent on only one 
trim condition simulation parameter. To control high bank 
turn maneuvers the simulation has been extended to automatically 
trim in a steady coordinated turn. The significance of this 
maneuver is the difficulty for the pilot of maintaining altitude 
under high g-loads. 

The simple linear pilot model represents the magnitude 
and phase response of a previously studied control-theoretic 
model of optimal human behavior. However, this is the most 
unreliable and variable part of the control loop. Suggestions 
for further research are given below. 

7.3 CONTROL DESIGN FOR THE ZOOM AND PUSHOVER 

Several conclusions can be made about design on a reduced 
order model. Significant feed forward terms result when there 


90 



a number of mid- frequency poles in which the designer is not 
strongly interested (e.g., the pilot model). Frequency shaped 
control methods are extremely helpful where the middle and high 
frequency model is obscured by model reduction. The feed forward 
term of the control into the sensor model is realized as a feed- 
back loop in the compensator, and can have a significant effect 
on compensator dynamics. 

Manipulation of models in a loop where each block is a 
linear system requires a flexible design program to form different 
block matrix equations which result from breaking the loop at 
different points. The abilities to quickly compute eigenvalues 
and frequency response plots are extremely helpful for catching 
modeling or equation errors, particularly once the physical 
interpretation of the model is obscured by model reduction. 

7.4 TECHNICAL PROBLEMS 

Since the validity of the linear models relies on a trimmed 
condition, the ability to get the SIMII simulation to trim near 
a desired flight condition can be a significant problem. Since 
the maneuvers of interest are all far from straight and level 
flight, a significant number of simulation runs were made to 
trim. A trimmed condition for the high bank turn necessary to 
obtain linear models for the level and windup turn maneuvers was 
not achieved during the contract. Hence, the linear control design 
methods have all been illustrated on the zoom and pushover maneuver. 
Changes to SIMII to help facilitate more efficient trim solution 
are described in the FTTC User's Guide (Appendix C). Further 
simulations, or identification not requiring trim, can provide 
these models. 

The controlled plant transfer function is approximately 

2 

.05/s . Yet increasing the pilot gain by 20. is too much to 
keep the low frequency dynamics from going unstable. Hence, 
additional work on pilot modeling and this gain adaptivity would 
be desirable. 
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7.5 RECOMMENDATIONS FOR FURTHER RESEARCH 


Although robust controllers can be designed with existing 
pilot models, their lack of fidelity does not parallel detailed 
knowledge of the aircraft and its control system. This situation 
can in some sense be remedied by changing the structure of the 
FTTC to feed forward directly into the autopilot as well as 
through the pilot by driving a display. The task can be separated 
by frequency, giving the pilot his most favorable frequency band 
and letting the autopilot handle low and high frequency control. 
Such a controller can be developed by frequency weighting these 
two portions of the gain inversely with the general shape of the 
pilot model in a wide enough range that one has effectively, 
"safely" inverted a poorly known pilot in the loop to get higher 
performance with the controller. 

Both the difficulty in achieving trim and all the engineering 
effort involved in modeling an MCS/CAS system can possibly be 
alleviated using linear identification and the known structure 
of the plant, identifying the aerodynamics involved in the 
maneuver and parameters of an assumed structure for the MCS/CAS. 

It would be very desirable to validate an FTTC for a 
difficult maneuver, such as a high bank turn in a manned simulator 
after sufficient optimization of the controller in the batch 
simulation. Data could be collected in such a study for identi- 
fication of the pilot model as well. 
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APPENDIX A 


FREQUENCY-SHAPING OF COST FUNCTIONALS 


Consider state-variable model of the form 

x = Fx + Gu + w , (A.l) 

with slate vector x, input vector u, and state-noise vector 
w. In steady-state LQG designs, a cost functional of the 
following form is minimized: 

r T 

J = lim if (x T Ax + u^Bu) dt . (A .2) 

T-*« JO 

A is i posit ive-semidef inite matrix and B is a positive- ■ 
definite matrix. In the frequency domain 

/ oo 

( x* Ax + u*Bu) dw , (A. 3) 

.00 

where x in Eq . A. 3 is the Fourier transform of time domain 
x. A anji B are the sarre as in Eq . A. 2. Note that Eq . A. 3 
speci' i es an equal penalty for states and inputs at all 
freque icies. 

1 le problem with the time-domain approach is that it is 
not pc ;sible to specify frequency-dependent penalties. There- 
fore, n Ref. [9], the concept of frequency-shaping was developed 
in the frequency domain. He can make A and B functions of 
freque cy in Eq. A. 3. A technique for designing control laws 
for ra _ Ional functional weightings has been developed [9], The 
procedure requires appending additional states in LQG design 
problems. Instead of showing the mathematical details of the 
design approach, we will demonstrate it by an example. 
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Example : Consider a one-mode, undamped oscillator with 

state equations 


*o-i: :] 


u 


(A. 4) 


We shall select a frequency-shaped weighting function of the 
form 


■/ 


(x*ax + u*bu) d , 


(A. 5) 


a = 


o> 2 + 4 


(A. 6) 


b = to + 1 . 


(A. 7) 


Defining 


z = - — — o => z + 2z = x , 
ju) + 2 

u = (ju)+l)u = u + u , 


( A . 8 ) 
(A. 9) 


the co ;t functional becomes 


0 

/ 


J = / ( z*z + u*u) dw 


= lim f J 

T-*.oo 1 J 


( z 2 + u 2 ) dt 


* (A. 10) 


The appended state equations are a combination of Eqs . A. 4, A. 8 
and A. 9: 


_d_ 

dt 


~ x~ 


1 — 
o 

o 

T~ 1 

0 

1 


x " 


" 0 " 

• 

X 

_ 

-4 0 0 1 


• 

X 

+ 

0 

u 


0 0-1 0 


u 


1 

z 


10 0 -2_ 


z 


_ 0 _ 


u 


(A. 11) 
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Standard LQG methods may be used to develop control laws for 
this modified problem. The designed controller will take the 
following form: 


u = C 1 x + C 2 x + CgZ + C 4 u . (A. 12) 

The implementation will use the equations obtained by substituting 
Eq. A. 12 into the last two terms of Eq . A. 11: 





(A. 13) 


Note that the controller is a low-pass filter. Thus, the region 
where the penalty on the state is small and the penalty on the 
control is high has reduced control activity. This behavior 
is typical of frequency-shaped control laws. 

Notes : 

1. The posit i ve-semidef inite matrices A and B, 
provided they are functions of to 2 , may always be 
decomposed to convert a frequency-shaping problem 
into a modified LQG problem. 

2. A standard Riccati equation needs to be solved to 
obtain the solution to the frequency-shaping 
problem. 

3. The resulting controllers are usually dynamic and 
have memory. In fact, frequency-shaping is achieved 
through this memory. 

4. The frequency shapes allowed by the techniques of 
Ref. [9] are a ratio of the polynomials in to 2 . A 
wide class of shapings may be approximated by such 
functions . 

5. The frequency-shaping techniques shown above for 
the control design problem may also be extended for 
state estimation (see the following section). 
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A.l ROBUSTNESS TO MODELING INACCURACY 


Effects of inaccurate open-loop models on the closed-loop 
performance is m maged by injecting minimum control power at 
natural frequencies of unmodeled modes. We will discuss proce- 
dures for controlling model inaccuracies in the high-frequency 
region, although similar techniques are applicable for other 
regimes. 

The high-frequency robustness can be accomplished by modi- 
fying the state of the control weighting. The state weighting 
A(joj) may be made a decreasing function of frequency. Three 
of the possible forms for A(jw) 

(1 ) A( jui) = -Tjj A , 

U) 

4 

w 

(2) A(ju>) - 2 2 2 ^ * 

(<jl) o + 0) ) 

. 2 . 2 .2 

( u- + u ) 

(3) A(ju)) = g 2 2 ^ ’ 

( + co ) 

To treat the frequency shaping of Eq . A. 14, we define additional 
states x as follows: 

x = x ( A . 17 ) 


are shown in the following: 


(A. 14) 


(A. 15) 


u>i > u) . 

1 o 


(A. 16) 


The performance index is 

.T 


1 f (-■ 

im = I (x 
■►00 


lim ^ / (x T lx + u^Bu) dt 

T+oo 


(A. 18) 


The control law will be of the form 


u = C^x + CgX 


(A. 19) 
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The control gain C 2 will ensure that the high-frequency 
response is minimized. This formulation resembles integral 
control. However, the weighting matrices used in the design 
are different . 

Spillover may also be reduced by placing high control 
weighting at high frequency. Examples of B(jw) which reduce 
high-frequency spillover are as follows: 


, 2 , 2 . 
(to 4a> 0 ) 


(1) B(jto) = B , 

to” 


(A. 20) 


(2) B(j u ) = 


2 . 2 
io^ + to 

2 . 2 B • W 1 < w o * 

10 + <0 


(A. 21) 


(3) B(jw) 


to 


[(CO 2 + %) 2 + 2ujw 2 ] 


B . 


(A. 22) 


The frequency-shaped parts of these weightings are reciprocals 
of the frequency-shaped parts of state weightings. To implement 
the weighting function of Eq . A. 20, we define a vector u as 


U + Ul U = lii u 

o o 


(A. 23) 


The performance index takes the form 

.T 


J = li 
T-*-® 


lm T J 

+-OD 


(x^Ax + u^Bu) dt , 


% 

{A. 24) 


The feedback control law will be 


U4u)U = tOU = U) (C.x 4 C 0 u) , 
O O 0 12 

or 


L 4 W o (l-C 2 )- = CjX . 


(A. 25) 
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This control law is shown schematically in Figure A.l. 

During maneuvers and in other missions, there is often a 
need to not excite a particular mode at w . Minimization of 
model errors at one frequency w o is achieved by: 

A(ju>) = (w 2 - w 2 ) 2 , (A. 26) 


or by 


B( j w) 



B 


2^2 


(A. 27) 



Figure A-l. A Frequency-Shaped Controller to Reduce 
High-Frequency Model Errors 

A. 2 ROBUST STATE ESTIMATION 

The discussion will again address the estimation of low- 
frequency states when the high-frequency modes are not modeled. 
Extensions to other cases is straightforward. * 

In the linear systems of Eq . A.l, a linear transformation „ 
can be used to place all errors cause by unmodeled modes into 
the measurement equations. Note that when high-frequency modes 
are not modeled in the state equations, measurement errors in 
the modified system are also at high freqency. Therefore, the 
lumped modeling error at high frequency represents a high- 
frequency noise measurement with many peaks. 
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To derive a filter with desirable high-frequency behavior, 
we consider optimal estimation as the output of an optimization 
problem with the following performance index: 


J 


1 im 

T-t-oc 


1 

T 



T 

w Ow 


v^Rv) 


dt 


(A. 28) 


Again, conversion to frequency domain gives the following 
performance index: 


J = 



(W*Q( jld)W 


— OC 


+ VRCjojJv) du . 


(A. 29) 


The problem of model invalidity at high frequency is solved by 
making Q and R functions of frequency. 


Since all errors associated with model truncation are 
incorporated in the measurements, we will consider frequency 
shaping only in R. High-frequency measurement spillover is 
controlled by increasing R(ju) at high frequency. We may 
select R(ju) as: 


R( ju>) 



+ u 2) 
2 


R . 


(A. 30) 


We can define a new measurement y as 


V = 


O' 


(A. 31) 


and 


— = Hx + v 

y 


where 


(A. 32) 


x + 


u x . 
o 


(A. 33) 
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The state estimator is of the form: 
x = Fx + Gu = K(y - Hx) , 
x = -lo q x + u> Q x + to q K ( y - Hx) , 

y = -« 0 y + oo Q y • (A. 34) 

This filter is shown schematically in Figure A-r2 . Note that 
a standard filter is obtained by setting K to zero and 
oj q /(s + a> o ) to one. In this formulation, filtered measurements 
are compared with filtered states. The eigenvalues of the closed- 
loop system are, in general, in the low-frequency region. Low- 
frequercy, closed-loop eigenvalues reduce excitation of high- 
frequency modes through feedback in the state estimator. The 
filter shapes may be modified to provide any desired roll-off. 

The fi ter, of course, becomes more complex as faster roll-off 
is obt -ined. Note also that the transfer function between the 
estima ed state and the measurements has at least two more poles 
than zros. This filter has wide applicability in the design 
of pra tical state estimators. 



Figure A-2 . A Frequency-Shaped Filter to Reduce 
High-Frequency Modeling Errors 
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A. 3 DISTURBANCE REJECTION IN FLIGHT VEHICLES 


Aircraft and rotorcraft are subject to high disturbance 
levels caused by engines or moving subsystems during in flight 
reconfiguration. In addition, turbulence, thrust variations, 
and density variations are external disturbance sources. These 
sources are nonwhite. Frequency-shaping methods may be applied 
to minimize the effect of these disturbances on the aircraft. 
Procedures for realizing disturbance rejection are shown in the 
fol lowing. 

Consider a configuration transition of the F-lll. The 
primary disturbance produced by the transition is at one frequency 
(typically the rotational frequency) and possibly its harmonics. 
Let the disturbance frequency be to. Suppose we wish to minimize 
the effect of this disturbance on output y. This disturbance 
rejection can be realized by including the following term in 
the performance index: 


.2 1 — 2^2 yT V • 


(A. 35) 


u) - CO ) 


Note that the output penalty goes to infinity at the disturbance 
frequency. Thus, the impact of disturbance on the output is 
minimized . 

The implementation of this frequency-shaped weighting 
requires definition of additional states as follows: 


x = 




(A. 36) 


The performance index will take the form 


J = 


lim 

T 


f T - 
T J (X 


- T A x + x T Ax + u T Bu ) dt , 


(A. 37) 
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and the control law is 


u = C^x + C 2 X 

The flowchart for this control law is shown in Figure A-3. 

The control formulation given above may be extended to 
include disturbance at other frequencies or over a frequency 
band. In each case, as shown in Figure A-3, there is memory 
in the control portion of the feedback. 

A. 4 SUMMARY 

The application of frequency-shaping methods to flight 
vehicles leads to a linear controller with memory. However, 
the additional states needed to emulate frequency-dependent 
weights increase controller order. Software needed for these 
controller designs is similar to that for standard LQG problems. 



Figure A-3. A Frequency-Shaped Controller to Eliminate the 
Effect of Disturbance at 00 on Output y 
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APPENDIX B 


MATRIX X : A DATA ANALYSIS, SYSTEM IDENTIFICATION CONTROL 

DESIGN AND SIMULATION PROGRAM 


Robert Walker, Charles Gregory, Jr., and Sunil Shah 


MATRIX X is an interactive software system to perform a com- 
plete cycle of steps starting from data analysis to system iden- 
tification, control design, simulation and evaluation. It is 
built on a user-friendly interpreter incorporating powerful 
matrix operations. State-of-the-art algorithms for linear 
system analysis, differential equation solutions, and Fourier 
transformation are included and can be called from the interpreter 
by simple commands. The package offers four major features 

(1) Powerful interpreter, simple command structure, good 
graphics capability with most "book-keeping" chores 
handled by the software, 

(2) State-of-the-art numerical algorithms, which allow 
solutions to high order problems, 

(3) User transparent file management with uniform and 
consistent format, and 

(4) Efficient implementation with a stack to require minimum 
memory and computation resources. 

Basic algorithms for linear system solutions, eigensystem 
decomposition (including reliable determination of the Jordan 
form), singular value decomposition (SVD), QZ decomposition, 
and matrix algebraic operations are implemented as language 
primitives. Many of the primitives are inherited from MATRIXES 
predecessor MATLAB, developed by Moler. 
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The interpreter can execute higher level constructs called 
macros and command files enabling construction of multi-level 
hierarchical structures. The command files are typically less 
than ten lines long and consist of interpreter primitives or 
other higher level constructs. Creation and modification of 
specific design and analysis procedures is easy because of the 
rich set of primitives and the hierarchical structure. 

A set of command files and macros provide specific environ- 
ments, such as data analysis and identification, control design 
and analysis, and simulation. 

Data analysis and identification can be performed very effi- 
ciently and easily in MATRIX^. Tied with a flexible graphics 
package, MATRIX X provides a production environment for batch and 
recursive identification methods. A universal interface with 
external simulations is provided to facilitate data transfer. 

Data can be censored, detrended and analyzed in MATRIX X> Batch 
procedures include the standard regression methods with analysis 
of variance and step-wise regression. State-space and nonlinear 
batch maximum likelihood procedures are also available. Recursive 
algorithms such as the recursive maximum likelihood, extended 
Kalman filter with Ljung's modification, and recursive instru- 
mental variable method are implemented. All covariance factori- 
zations and updates are in U-D form for numerical reliability. 
Non-parametric batch and semi-batch methods using the FFT are 
provided for auto/cross covariances/spectras. Lattice updates 
for ARMAX systems are available. Adaptive control algorithms 
for multivariable systems using U-D updates can be designed 
using simple commands. 

In MATRIX X , control design can be based on any of the 
following: 

(a) Linear Quadratic Gassian (LQG) approach, 

(b) Methods based on A-B invariant subspaces, 
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(c) Eigenst ructure assignment and zero placement, and 

(d) Classical methods such as Nyquist and Bode. 

For the LQG problem, the algebraic Riccati equation is 
solved from extended Hamilton equations avoiding inverses, which 
are troublesome in the singular case. The equations are row 
compressed with an orthogonal transformation followed by the 
QZ pencil decomposition and a backward stable ordering of the 
eigenvalues . 

Meaningful extensions to LQG methods require inclusion of 
dynamics of reference, disturbances, sensors, and actuators. 
Appending of dynamics in frequency-shaped control design or 
model-following techniques involves forming augmented equations, 
which is easily accomplished with MATRIX X primitives. Use of 
frequency-shaped costs, with singular value plots for robust- 
ness evaluation, allow incorporation of engineering judgment 
in the control design. 

Evaluation tools for linear systems include frequency 
response, power spectral density plots, time reponses, trans- 
mission zeroes and individual transfer function zeroes. The 
principal vector algorithm (PVA) primitive for numerically 
reliable extraction of the Jordan Form (with discriminatory 
rank deflation of root clusters) is very useful in modal analy- 
sis of open-loop systems of vehicles and structures. PVA 
permits computation of residues or partial fraction expansions 
of multivariable systems. Extraction of the Kronecker indices, 
supremal (A,B) invariant subspaces and the Kalman decomposition 
is performed using the SVD, QR and the QZ algorithm. 

Simulation of stiff differential equations is performed 
with Gear's software. Linear system responses are computed 
efficiently from the residues. 

Use of "chopped arithmetic," i.e., using various effective 
machine word-lengths can provide performance evaluation of 
on-board small word-length control system implementations. 
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Numerical reliability and stability are important in all of 
the analysis environments described above. Primitive matrix 
operations are based on the best available numerical software 
drawn from EISPACK, LINPACK, and recent research in numerical 
techniques . 

Reporting and control of numerical errors is performed in 
several ways: Estimates of the problem condition number and 

the algorithms condition number are available for many of the 
primitives in MATRIX X . When requested, MATRIX X provides an 
estimate of the problem condition, the algorithm condition and 
the solution accuracy when such an analysis has been presented 
in the current literature. If an error analysis has not been 
performed rigorously, condition estimates using varying word 
length arithmetic and perturbed initial data can be requested. 

The program is coded in ANSI-77 FORTRAN and will run on 
any system allowing interactive execution of FORTRAN. The few 
machine dependent features are implemented separately for most 
popular operating systems. An implementation with 250K double- 
precision complex elements is available for VAX 11/780. System 
orders up to 150 can be handled easily within the storage 
capacity of the computer and numerical accuracy of the algorithms. 
Most algorithms access memory in a linear fashion, so that page- 
swapping is not unduly increased in the virtual memory implemen- 
tat ion . 

In summary, MATRIX X provides a system which minimizes engi- 
neering and programming resources required for the complete 
cycle of system identification, control design, and validation, 
yet provides the flexibility needed in an applied research 
environment . 
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APPENDIX C 


FLIGHT TEST TRAJECTORY CONTROLLER USERS* GUIDE 


AUTOTRIM CAPABILITIES 

In the version of the NASA-Dryden Simulation SIMII modified 
for flight test trajectory control (FTTC), the autotrim function 
is quite important since the linear models depend on achieving 
trim. The previous four autotrim options are available with two 
new autotrims for steady coordinated turns. Level turns require 
a type of dynamic trim, since there are constant body angular 
rates, rather than zero angular rates. One of the turn trims 
specifies attitude, effective roll angle and angle of attack, 
adjusting the turn rate and velocity to achieve trim. The second 
turn trim is vice versa, specifying turn rate and velocity, 
adjusting the effective bank angle and angle of attack to achieve 
trim. A summary of constrained and adjusted variables for all 
six trim conditions is shown in Table 1 below. 


Table 1 SIMII Trim Options 


Trim Type 

Program 

Mnemonic 

Constrained 

Variables 

Adjusted 

Variables 

Normal 

NORMAL, also 
'febbbfeb ' 

V 

PLA,a,B 

Throttle 

THRP 

PLA 

V,a,B 

Angle-of-Attack 

ALP 

a 

PLA , V , 6 

Sideslip 

BTA 

6 

a ,V,4> 

Turn rate 

PS I DOT 

*o’ V 

PLA , a , 
4> L ( = >^ ,6 ,<f>) 

Effective roll 
angle (angle of 
lift vector from 
the vertical) 

PHILFT 

a, <I> L (=><!>, 6 ,<)>) 

PLA , V , 4 1 
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The gains which adjust the unconstrained variables at each 
autotrim iteration can be entered as input. Table 2 shows the 
appropriate user input variable and the associated default value 
for zero input. 


Table 2 Autotrim Adjustment Gains and Input Parameters 



UX( 79 ) 

Starting value for PLA (Default = 0) 

UX(80 ) 

Maximum number of iterations (10,000) 

UX(89) 

0 

ip o (Turn rate trim option )=> p Q , q Q , r Q 

UX( 90 ) 

<j>k (Effective bank angle option), 

a => il> ,8 , d> 
o r o’ o ,Y o 


The trim values can be preset to desired values without 
going through the autotrim iterations with the subroutine TRMSET. 
Values of the UX array are used as input. The association of 
variables and locations of UX are documented in this subroutine, 
included in Appendix D. 

A summary of the user logical flags useful in FTTC design 
and evaluation are given below in Table 3. 
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Table 3 SIMII-FTTC Option Flags 


Logical Input Flag 

Option 

LX (20) = T 

Find Stability Derivatives (see 
Aerodynamic Linearization section) 

LX (21) = T 

Stop SIMII after finding stability 
derivatives 

LX (22) = T 

Bypass trim but set all trim values 
(invokes subroutine TRMSET) 

LX ( 23 ) = T 

Print out intermediate autotrim 
convergence data 

LX ( 24 ) = T 

FTTC closed loop simulation (additional 
inputs required (see Tables 4 and 5) 

LX( 25) - LX( 28 ) = T 

Invoke a pilot model for, respectively, 
the Aileron, Elevator, Rudder, and 
Thrust commands 


AERODYNAMIC LINEARIZATION 

When the simulation SIMII has been trimmed, all the stability 
derivatives can be requested as follows: 

LX(20) = T Computes all the stability derivatives, 

LX(21) = T Stops the simulation after the stability 

calculations 

The six aerodynamic force and moment coefficients and the inde- 
pendent variables that affect them are shown in Table 5-2, the 
zoom and pushover trimmed example. Both stability axis and 
body axis variables are shown. The equations in SIMII are imple- 
mented with the force coefficients in the stability axes and 
moment coefficients in body axes. The linearization procedure 
perturbs the appropriate combinations of variables in stability 
axes to get the body axes derivatives. The effect of a center- 
of-gravity shift is included and the divided differences of the 
nondimensional coefficients are converted to dimensional stability 
derivatives before output. 


Ill 












The subroutine STABDR outputs the trimmed nondimensional 
coefficients along with the dimensionalizing coefficients (for 
checking total force and moment balances), sets up consistent 
velocity perturbations, and calls the routine SDCALC for each 
state variable of interest. SDCALC saves the initial coefficients 
in another six element vector, perturbs the passed state variable 
of interest, recalls the SIMII routine CCALC to find the per- 
turbed six coefficients, finds the divided differences, appro- 
priately dimensionalizes these stability derivative approximations, 
and prints out the state variable with its six stability deriva- 
t ives . 

FTTC INPUT AND OUTPUT 

FTTC simulation is invoked by inputting LX(24) = T, with 
the necessary inputs summarized in Tables 4 through 6. Of the 
various blocks of the flight test trajectory controller (see 
Figure 6-1), the FTTCOM common block is used for linear models 
and the UX array is used to store the outputs of the various 
stages of the controller. Thus, time histories of the internal 
variables of the controller can be output for plotting. 

Table 4 FTTC Input and Output Variables 


Variable 

Storage Location 

Aircraft Model Trim Values 

US(ll) - US(25) 

Sensor Noise Standard Deviations 

US( 26 ) - US(40 ) 

Trajectory Constants 

US(41 ) - US( 55) 

Compensator to SIMII Unit Conversions 

US( 56 ) - US( 59 ) 

Sensor Outputs 

UX(ll) - UX ( 2 5 ) 

Trajectory Command Generator Outputs 

UX ( 26 ) - UX (40) 

Compensator Inputs 

UX (41) - UX( 55 ) 

Compensator Outputs 

UX( 56 ) - UX ( 59 ) 
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The input for the compensator and pilot models is accom- 
plished by FTTCIN, a subroutine called where the control input 
data is read, CNTRIN. Thus the FTTC input comes immediately 
after Group 11 data cards, the autotrim variables. The necessary 
FTTCIN cards are given in Table 5 below. 


Table 5 FTTC Linear Model Input 


Card 

Column 

Format 

Description 

1 

1-5 

6-10 

215 

Dimension of the compensator dynamics 
matrix (F c ), and the number of com- 
pensator inputs (i.e. , NZ the number 
of sensors) 

2 

1-60 

3D20. 13 

F c , G c , H c , J c , by rows (G c by 
columns), a new card for each row 
(column). Note that there are four 
compensator outputs for the four 
principal controls 

3 

1-5 

15 

Dimension of the first pilot model 
(e.g., Fp A ). 

4 

1-60 

3D20. 13 

Fp a , GpA, Hpa, JPa by rows (G by 
columns), a new card for each row 
(column). Note that each pilot model 
is an independent SISO system 

5 

• 

• • • • 

• 

Additional pilot models 

j 




Immediately following the FTTC linear system input, the 
necessary FTTC constants are input as shown in Table 6. 
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Table 6 FTTC Input Parameters 


Card 

Type 

Column 

Number 

Format 

Description 

Variable 

Locations 

1 

1-5 

6-10 

215 

Sensor and pilot model 
sampling intervals (i.e., 
multiple of simulation 
integration step) 

NSAMPZ , 
NSAMPP 

2 

1-50 

5E10 . 3 

Aircraft trim values, 
y T (NZxl) trim values plus 
four control trim values 

US(ll) - 
US(25) 

3 

1-50 

5E10 . 3 

Sensor noise standard devia- 
tions (assuming uncorrelated 
sensors) 

US( 26 ) - 
US ( 40 ) 

4 

1-5 

15 

The number of trajectory 
parameters 

NTRAJP 

5 

1-50 

5E10.3 

Trajectory parameter 
constants 

US( 41 ) - 
US(55) 

6 

1-50 


Compensator to SIMII control 
unit conversions 

US( 56) - 
US(59) 
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APPENDIX D 


FORTRAN LISTINGS FROM SIMII-FTTC 

The new or significantly changed subroutines of the NASA 
DRC F-15 SIMII simulation are listed by function below. Their 
calling sequence and function are shown in Figure 6-2. The 
FORTRAN listings follow chronologically. 

Autotrim Capabilities 

AUTOTR Adjust variables to achieve desired initial 
condition 

TRMCHK Check force and moment balances 

TRMSET^ Set all trim values to desired variables 

Aerodynamic Linearization 

STABDR 1 Perturb each state variable of interest 

SDCALC^ Compute six stability derivatives for each 
state variable 

FTTC Simulation 

FTTC^ Integrate flight test trajectory control law 

SENSOR-^ Extract SIMII variables and add sensor noise 

TRAJGN 1 Trajectory command generator 

PROPAG Propagate discrete compensator and pilot models 

MATRIXx Control Design Program Interface 

MTOSIM Translate MATRIXx data files with compensator 

and pilot linear models to SIMII-FTTC compatible 
file format 

SIMTOM Translate SIMII time history file to MATRIX^ 
compatible file format for plotting 

^Listing attached to this report 
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SUBROUTINE TRMSET 

IMPLICIT DOUBLE PRECISION (A-H»0-Z> 

LOGICAL LX 

C TRMSET SIMULATES A TRIMMED CONDITION BY SETTING THE VALUES 

C FROM UX<21)-UX<63> 

C 

COMMON /VARUAT/ UX< 90 ) » LX ( 50 > 

COMMON /DRVOUT/ F(13)»DF(13> 

COMMON /SIMQUT/ AMCH * QBAR * GMA * BEL * UB > VB * t*B * VEAS » VC AS 
COMMON /CONTRL/ D A . BAR * DAL > DFF * DH> DHR , DHL * DT * DR * DSB * THRST * THR , OHGE 
COMMON /CLCOUT/ CL*CM,CN*CD*CLFTtCY 
COMMON /C0NP03/ Fl(9)*DSBF*DFP*r2{4> 

COMMON /ENGDSF7 POSNZL , POSNZR » RPHL » RF'MR » TEMPL iTEMPR , 

. FFLOHL»FFLOHR»FQNTYi 

LOGICAL FOUT » NOWC 

COMMON /ENGINE/ FONT V » FFLOW » D »FOUT > NOUC » 3INET A » CQSETh 
DATA DTR/.0174532/ 


UX 21-46 ARE F.DF 

DO 10 I=lf 13 
Fa) - UX ( 20+1 ) 

10 DFil) = UXC33+I > 

CONVERT TO INTERNAL UNITS ON ANGLE STATES 

DO 20 1-6*10 
20 F ( I ) - FiUlDTR 

CONPOS VARIABLES UX 48-56*57-60 

DO 30 1-1*9 
30 F1(I) - UX047+I) 

DO 40 1-1*4 
40 F2(I) = UX(5o+I) 

CONTROL POSITIONS UX oi-64 

DA - UX(61)tDTR 
DH = UX(62)tDTR 
DT = UX(63)tDTR 
DR - UX(64XDTR 

ENGINE PARAMETERS UX 65-68 

FFLOU UXvo5) 

THRST -- UX(66) 

FFLOUL - UX( 67 ) 

FFLOMR = UX(68) 

AMCH ^ UX ( 7 i 7 
VEAS ^ UX (72) 

RETURN 

END 



SUBROUTINE STABDR 

IMPLICIT DOUBLE PRECISION (A-HfQ-Z) 

LOGICAL NORMAL , LONGI T > LATDIR f LEVEL 

CONNOR /SIMOUT/ AMCH f QBAR f BNA , DEL f UB r VB f WB f VEA3 f VCAS 

COMMON /ALTFUN/ AfRHOfGfPA 

COMMON /SIHTYF'/ NORHAL f LONGI T > LATDIR > LEVEL 

COMNON /HRVOUT/ F(13) »HF(13) 

COMMON /CLCOUT/ CLfCMfCNfCDfCLFTfCY 

COMMON /CQNTRL/ DA f DAR f DAL f DFF f DH f DHR f DHL f DT f DR f DSB > THRST » THR » OMGE 
COMMON /SIMACC/ AXfAYfAZfANXfANYfANZfAN 

COMMON /FMCOEF/ FM(39) 

COMMON / TRIGFN/ S INALP f COSALP f SINBTA f COSBTAf SINPHI f COSPHI f SINPSI f 

. COSPSIfSINTHAfCOSTHA 

COMMON /DATAIN/ SfBfCBARfAHSSfAIXfAIYfAIZfAIXZfAIXE 

COMMON /3TABCN/ QSfQSB 

EQUIVALENCE (T f F( 1))f 

(F f F! 2) ) f (Q f F< 3) ) f (R f F( 4>>f 

. (V f F( 5 ) ) f ( ALP f F( 6))f(BTA f F( 7>>f 

. (THA f F( 8 ) ) f < PSI f F( 9) ) f (PHI , F(10))f 

. iH f F(il)MX f F(12))f(Y f F(13>) 

STAB ILI TY VARIABLE PERTURBATION VALUES 

DATA CVDEL f ALPDEL f BT ADEL f PDEL f QDEL f RBELf UDEL f VDEL f HDEL f ANZDEL f 
. ADOTDL f DSBDEL f DADEL f DRDEL f DHDEL f BTDEL/16* .001/ 

DATA C VNAM f AN AM f BNAM fF’NAMf QNAM f RN AH f UNAM f VNAM f HNAN f AZNAM f ADNAM f 
. DSBNAMf DANAM , DRNAH » DHNAM f DTNAH/ ’ VELO ' f ' ALPA ' f ' BETA VP 'f 

. 'Q 'f R 'f'U VV 'f'« ‘f'ANZ VALFDVDSB 'f 

. 'DA VDR 'f'DH VDT 7 

FIND DIMENSIONAL CONSTANTS 

OS ^ OBARtS 


□SB ^ Q3*B 

WRITE ( 3 » 8900) AHSSf QBAR f QSfQSB 
3900 FORMATOV AMSSfQBARfQSfQSB = 'f4(1Xf1PE12.5)f/) 

WRITE\3f891Q) CLfCMfCNfCDfCLFTfCY 

8910 FORMAT ( ' NONDIMENSIONAL CLfCMfCNfCDfCLFTfCY = 'f/6(1Xf1PE12.5)f/> 
UR I TE < 3 f 8920 > <FM< I / f 1=1 f39) 

3920 FORNAT(/V CONSTITUENT COEFFICIENTS-SEE COMMON FHCOEF'f/f 
. (6(1Xf1PE12.5))) 

C SET UP VELOCITY PERTURBATIONS 

CVDEL = VtVDEL 
UDEL ^ CVDEL 
VDEL = CVDEL 
UDEL = CVDEL 
AMDEL = CVDEL/A 


C NOTE THAT STABILITY DERIV. ARE IN THE APPROPRIATE AXES 

URITE ( 3 f 9000 ) 

9000 FORMAT ( ' 1DIHENSI0NAL STABILITY DERIVATIVES FOR' f/f 

. 9Xf 'CLB0DY'f7Xf 'CM'fIIXf 'CNBODY' f7Xf'CD'f11Xf'CLFT'f9Xf ' CYBOBY ' ) 

C FIND LONGITUDINAL DERIVATIVES 

IF(LATDIR) 60 TO 100 
C - — ALPHA DER. 

CALL SDCALC ( CL f ALPDEL f ALP f ALPDELf VI f D1 f , FALSE , f 0 f ANAM ) 



c 


----- C HER. 

CALL SDCALC ( CL > QDEL » Q , DDEL »V1»D1» .FALSE . »0» QNAM ) 

U HER. 

CALL SDCALC i CL » MDEL »V»Dl»ALP»D2r .TRUE. »2» WNAH) 

C — - ANZ HER. 

CALL SDCALC < CL » ANZDEL i ANZ i ANZDEL »V2»D2» . FALSE . 1 0 , AZNAH > 

C ALPHADOT DER. 

CALL SDCALC ( CL > ADQTDL » DF ( 6 ) , ADOTDL » V2 >D2» . FALSE , , 0 > ADNAM ) 

C DERIVATIVES FOR BOTH LONGITUDINAL AND LATERAL DER. 

100 CONTINUE 
C VELOCITY DER. 

CALL SDCALC ( CL » CVDEL i V » CVDEL » AHCH» ANDEL » . TRUE .,0> CVNAH ) 

C U DER. 

C ALL SDCALC ( CL ,UDEL » V , D1 » ALP » D2 » . TRUE . » 1 » UNAH ) 

C - — DH DER. 

CALL SDCALC ( CL . DHDEL . DH » DHDEL » DSB » DHDEL » . TRUE . » 0 » DHNAM ) 

C COHPUTE STRICTLY LATERAL DERIVATIVES 

IF iLONGIT) GO TO 200 
C BETA DER. 

CALL SDCALC ( CL .BTADEL , BTAiBTADEL » V2 » D2» . FALSE . > 0 . BNAri ) 

C — - P HER. 

CALL SDCALC ( CL ,PDEL » P » PDEL i V2 »D2 , .FALSE . rO.PNAH ) 

C R DER. 

CALL SDCALC ( CL > RDEL > R > RDEL . V2 > D2 t .FALSE. »OiRNAH) 

C V DER. 

CALL SDCALC ( CL » VDEL t BT A . D1 » V2 » D2 . . FALSE . » 3 , VNAM ) 

C — - DA DER. 

CALL SDCALC ( CL » DADEL > DA » DADEL >V2»B2>. FALSE . r 0 » DANAH > 
t DR DER. 

CALL SDCALC ( CL » DRDEL » DR . DRDEL > V2 » D2 . . FALSE . » CrBRNAH ) 

C - — DT DER. 

CALL SDCALC ( CL r DTDEL > DT » DTDEL » V2»D2» . FALSE , . 0 1 DTNAH) 

200 RETURN 

END 
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c 

SUBROUTINE 30CALC ( C» DELTA* VAR1 r DELT A1 » VAR2 f DELTA2 * HULVAR » IBODY 
. »NAME) 

IMPLICIT DOUBLE PRECISION <A-H*Q-Z) 

LOGICAL HULVAR 
INTEGER HAHEj IBODY 

C — HULVAR - F => ONLY ONE VARIABLE IS PERTURBED 
C T => PERTURB TWO VARIABLES 

C 

C IBODY = 0 => USE STAB. AXES FOR STAB. VAR. AND BODY AXES FOR BODY VAR. 

C 1 => PERTURB VT AND ALPHA FOR U 

C 2 => PERTURB VT AND ALPHA FOR U 

C 3 => PERTURB VT AND BETA FOR V 

C 

DIMENSION C i 6 ) » CSAVE <6>»DC(6) 

COHHON /CGSHFT/ DELX r DELY » DELZ 

COHHON /DRVOUT/ K(13)*DFU3) 

COHHON /DATAIN/' 3*B*CBAK*AHSS*AIX*AIY»AIZ*AIXZ,AIXE 

COHHON /STABCN/ QS»QSB 

COHHON /TRIGFN/ SINALP*COSALP»SINBTA*COSBTA»SINPHI,COSPHI*SINPSI* 

COSF’SI rSINTHAiCOSTHA 
EQUIVALENCE (V »F<5)) 


C REINITALIZE AND SAVE COEFFICIENTS 

CALL CCALC 
DO 10 1-1 » 6 
10 CSAVE (3) -- CCD 

C PERTURB VAR TABLES (ONE OR TWO) 


IF( IBODY .EG, 0) GO TO 15 
GO T0( 11 *12*13) *IBODY 

C PERTURBATION FOR U — NOTE THAT V IS TOTAL V(VT) 

11 DELTA1 = COSALPIDELTA 

DELTA2 - -3INALPIDELTA/ ( VtCOSBTA) 

GO TO 15 

12 DELTA1 = SINALPIDELTA 

DELTA2 = COSALPIDELTA/ (VtCOSBTA) 

GO TO 15 

13 BELTA1 = DELTA/COSBTA 
15 VAR1 = VAR1 + DELTA1 

IF (.NOT. HULVAR) GO TO 20 
VAR2 - VAR2 + DELTA2 
20 CALL CCALC 

C FIND DIVIDED DIFFERENCES 

DO 30 1=1*6 

30 DC ( I ) = ( C ( I ) -CSAVE ( I ) ) /DELTA 

C PREPARE TO ADD EFFECT OF C. G, SHIFT 

TERH = ( DC ( 4 ) tCGSBT AISI NALP + DC(6)*SINBTA*SINALP + DC(5)*C0SALP)I 
. 03 

TERH3 = QSI ( DELXI ( DC ( 4 ) ISINBT A - DC(6)»C0SBTA) - 
. DEL Y* ( DC ( 4 ) ICOSBT AICOSALP t DC(6)*SINBTA*C0SALP - DC(5)*SINALP)) 

C HAKE DERIVATIVES DIMENSIONAL 

DC(1) = (DC(l)IQSB I TERHIDELYJ/AIX 
DC(2! (DC(2)*QSICBAR - TERH*DELX)/AIY 
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DCi3) - (DC(3)tflSB + TERH3)/AIZ 
DC( 4) = -DC ( 4 ) KQS/AHSS 
DC ( 5 > = -DC(5)*QS/AHSS 
BC<6> = DC ( 6 ) tQS/AKSS 

C TRAMSFORM CL AMD CD TO BODY AXES 

IFIIBODY .EO. 0) 60 TO 40 

DCL8 = COSALPtDC ( 5 ) - SINALP»DC(4) 

DCD = SINALPIDC(5) 4 C0SALP*DC(4) 

DC ( 5 > = DCLB 
DC (4) = DCD 

C RESTORE THE INITIAL VALUES 

40 DO 50 1=1 f 6 
50 C(I) = CSAVE(I) 

VAR1 = VAR1 - DELTA1 
IF(MULVAR) VAR2 = VAR2 - DELTA2 

C PRINT OUT STABILITY DERIVATIVES 

URITE(3»9000) NAME i DC 
9000 F0RMAT(lX»A4i6(lX»E12«5>) 

RETURN 

END 



SUBROUTINE FTTC(ITIHE) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

INTEGER ITIHE 

COMMON /FTTCOM/FTT ( 200 ) * IFCI 3 ) » IGC( 3 ) * IHC( 3) , I JC ( 3 ) * 

, IFPA(3)*IFPE(3MFPR(3)*IFPT(3bIWK(3)* 

. NSAMPZ fNSAMPP 

*•» 

L 

C FTTC IMPLEMENTS A FLIGHT TEST TRAJECTORY CONTROL LAW AT 
C INTEGRATION INTERVAL I TIME 

t 

C COMPENSATOR AND PILOT MODELS ARE STORED IN ARRAY FTT WITH 
C POINTER INDICES 

C I.E. IFCd) ROW DIMENSION OF FC (COMPENSATOR DYNAMICS) 

C IFC<2> COL DIMENSION OF FC 

C IFC(3) LOCATION IN FTT WHERE FC IS STORED 

C 

c 

c 

c IFF! ( 1 5 ROW DIM. OF FPTHtPILQT DYNAMIC MODEL FOR THRUST) 

C IFPK2) COL DIM. OF FPTH 

C IFPK3) LOCATION IN FCC WHERE FPTH BEGINS 

C NOTE! FOR THE FIVE LINEAR MODELS THE CONVENTION IS TO STORE 
C F»G»H»J»X»Y (MODEL STATE AND OUTPUT) IN COMMON FTT 

C 

C NSAMPZ TIME INCREMENT FOR UPDATING MEASUREMENTS 

C N3AMPP TIME INCREMENT FOR UPDATING PILOT MODEL 

C 

C AIRCRAFT MODEL TRIM VALUES 

C SENSOR NOISE STANDARD DEVIATIONS 

C TRAJECTORY CONSTANTS 

C UNIT CONVERSIONS (COMPENS. TO SIHII) 

C 

C SENSOR OUTPUTS 

C TRAJECTORY COMMAND GENERATOR OUTPUTS 

C COMPENSATOR INFUTS 

C COMPENSATOR OUTPUTS 

C 

COMMON /CONPOS/ DAP f FLAT » DATR IH i DEP r FLON » DETRI M » DRP » FPED » DRTR IM r 
. DSBP»DFP>FLAPLfPLAPRf THSL r THSR 

COMMON /VARDAT/ UX(90)*LX(55) 

COMMON /VARSIM/ US(90)fLS(50) 

C 

C INITIALIZE INPUT/OUTPUT HISTORY VARIABLES 

IFdTIHE ,GT. 1) GO TO 10 
DO 5 1=11*59 
0 UX(I) = O.DO 
10 CONTINUE 


US(11)-US(25) 

US(26)-US(40) 

US(41)-US(55) 

US(5o)-US(59) 

UX(11)-UX<25) 

UX(26)-UX(40) 

UX(41)-UX(55) 

UX(56)-UX(59) 


C GENERATE TRAJECTORY COMMANDS 

C 

IF(H0D(ITIME-1 jNSAMPZ) ,EQ. 0) CALL TRAJGN(US(41)*UX<26>) 

C UPDATE MEA3 

C 

IF ( MOD (I T I ME- 1 1 NSAMPZ ) ,EQ, 0) CALL SENSOR ( UX ( 1 1) >US( 11 > *UX(26) * 
. UX(41)»US(26)*NZ) 

C 

122 C PROPAGATE THE COMPENSATOR 



L 


30 IX=IJC(3> + IJC(l)tIJC(2) 

I Y=IX + IFC(l) 

CALL PR0FA6 ( FTT ( I FC ( 3 ) ) fIFC( 1 ) fFTT(IGC<3>) »IGC(2) fFTT(IHC(3) ) * 

. IHC( 1 ) »FTT < I JC<3> > ? FTT (IX)»UX(41)»FTT(IY) »FTT ( 2UK<3) > ) 

C 

C SAMPLE AMD HOLD THE COMPENSATOR OUTPUT 

C 

IF< MOD< ITIME-lf NSAMPP) ,NE. 0) GO TO 40 
DO 35 1=1 »4 

35 UX(55+I > = FTTdY-l+I) 

C 

c FROPAGATE THE PILOT HODELS 

C 

40 IF(,NOT. LX<25> ) GO TO 50 

C AILERON PILOT MODEL 

It = IFPA(3)+IFPA(1)*IFPA(2) 

IH = IG + IFPA(l) 

IX = IH + IFPA(l) + 1 
IY = IX + IFPA(l) 

CALL F ROPAG ( FTT (I FPA ( 3 ) ) » IFPA (l)»FTT(IG)tl»FTT<IH)»l»FTT(IX-l)» 
. FTT(IX)fUX(56)fFTT(IY)fFTT(IHK(3))) 

C ADD A/C TRIM TO COMP. OUTPUT 

DAP = FTT(IY)*US<56) + US(lltNZ) 

C 

50 IF ( .NOT . LX<26)> GO TO 60 

c ELEVATOR PILOT MODEL 

IG = IFPE(3)+IFPE(l)tIFFE(2> 

IH = IG t IFPE(l) 

IX = IH + IFPE(l) + 1 
1Y = IX + IFPE(l) 

CALL PRQPAG(FTT(IFPE(3))fIFPE(1>fFTT<IG)f1fFTT<IH)f1fFTT(IX-1)f 


. FTT(IX)iUX(57)»FTT(IY)fFTT(IUK(3))) 

C ADD A/C TRIM TO COMP. OUTPUT 

DEP = FTT(IY)*US(57) + US(12+NZ) 

C 

60 IF ( .NOT . LX<27) ) 60 TO 70 
C RUDDER PILOT MODEL 


IG = IFPR(3)+IFPR(1)$IFPR(2) 

IH = IG + IFPR(l) 

IX - IH + IFPR(l) + 1 
IY - IX t IFPE(l) 

CALL PRQPAG(FTT(IFPR(3))fIFPR<1)fFTT(I6)f1fFTT(IH)f1fFTT<IX-1)f 


. FTT (IX) fUX(58) fFTT(IY) fFTT(IHK(3) ) ) 

C ADD A/C TRIM TO COMP. OUTPUT 

DRP = FTT(IY)*US(58) + USU3+NZ) 

C 

70 IF ( .NOT * LX (28) ) GO TO 80 
C THRUST PILOT MODEL 


IG = IFPT(3)+IFPT(1)*IFPT(2> 

IH - IG P IFPT(l) 

IX - IH + IFPT(l) P 1 
IY - IX P IFPT(l) 

CALL PR0PAG<FTT(IFPT<3))fIFPT(1)fFTT(IG)f1fFTT(IH)f1fFTT(IX-1)f 
. FTT(IX)fUX(59)fFTT(IY)fFTT(IHK(3))) 

C ADD A/C TRIM TO COMP. OUTPUT 

THRF = FTT(IY)*US(5?) P US(14PNZ) 
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80 CONTINUE 


C 

RETURN 

C 

ENTRY FTTCIN 
C 

C READ COMPENSATOR MATRICES 

C 

READ(lriCl) IFC(1)»IGC<2) 

101 FORMAT! 1615) 

IFC<2) ^ IFC(l) 

IGC(l) = IFC(l) 

IFC(3) = 1 

IOC (3) = 1 + IFC!1)«FC<2) 

IHC( i > = 4 
I HC ( 2 > = IFCU) 

IHC(3) = I GC ( 3 ) +IGC( 1 ) *IGC ( 2 ) 

IJC(1> - A 
I JC ( 2 ) = IGC(2) 

IX<3) = IHC<3) I IHC(1)*IHC(2) 

NA -- IFC(l) 

NZ = I GC ( 2 ) 

CALL RE ABM A ( FTT f N A f NA f NA f » TRUE . f ' ( 3D20 . 1 3 ) ') 

CALL READMA (FTT ( IGC ( 3 ) ) f NAf NA f IGC ! 2 ) f .FALSE . f ' ( 3D20 , 13 ) ') 

CALL READMA ( FTT ( IHC ( 3) ) f 4 f 4 f NAf . TRUE . f ' ( 3D20 .13) ') 

CALL READMA(FTT(IJC(3))f4f4fIGC(2)f»TRUE.f '(3D20.13) ') 

C 

WRITE! 3» 103) 

103 FORMAT!/ i ' FTTC COMPENSATOR MODEL') 

CALL WRITNA!FTTfNAfNAfNAf 'FC ') 

CALL WRITMA!FTT<IGC<3))fNAfNAfIGC(2)f'GC ') 

CALL WRITMA!FTT(IHC(3>)f4f4»NAf'HC ') 

CALL WRITMA!FTT(IJC!3))f4f4fIGC!2)f 'JC ') 

C 

WRITE(3f8900) (LX(24+I)fI=1f4) 

3900 FORMAT !/f' PILOT MODELS ARE ACTIVE FOR THE FOLLOWING CONTROLS- 
•SXt'DAP ' fL1f/f5Xf'DEP ' fL1f/f5Xf'DRP 'fL1f/f 
i5Xf 'THRP 'fLD 

1ST = IJC!3)FIJC(l)*IJC(2)+IFC(l)+4 

C 

IF!, NOT, LX(25) ) 60 TO 90 
WRITE(3f9010) 

9010 FORMAT'/f' AILERON CONTROL PILOT HODEL MATRICES') 

READ! If 101 ) IFPA(l) 

NA = IFPA(l) 

IFPA12) = NA 
IFPA<3) = 1ST 

CALL READMA (FTT! 1ST) fNAfNAfNAf, TRUE, f '(3D20, 13) ') 

CALL WRITMA! FTT! 1ST) fNAfNAfNAf 'FPA ') 

N2 = NAINA 

CALL READMA(FTT<ISTtN2)»NAfNAf1f, FALSE, f'(3D20. 13) ') 

CALL WRITMA! FTT (ISTFN2)fNAfNAf1f'GPA ') 

CALL RE ADMA ( FTT ( I ST+N2+NA ) f 1 f 1 f NA f , TRUE , f ' ( 3D20 .13) ') 

CALL WRITMA(FTT(IST+N2+NA>f1f1fNAf'HPA ') 

CALL READMA ( FTT (1ST +N2+2INA ) f l f l f 1 f ♦ TRUE , f ' ( 3D20 .13) ') 

CALL WRITMA<FTT(IST+N2+2*NA>f1f1f1f'JPA ') 



1ST » IST+N2+3INA+2 


C 

90 IF(.N0T. LX<26> > GO TO 100 
URITE(3»9020) 

9020 FORMAT (/j ' ELEVATOR CONTROL PILOT MODEL MATRICES’) 

READ! If 101/ NA 
IFPE(l) x NA 
IFPE ( 2 > = NA 
IFPEv 3) = 1ST 

CALL READMA(FTT(IST)»NA»NA»NA» • TRUE . » ' ( 3D2C ♦ 1 3 > ') 

CALL URITHA(FTT ( 1ST) >NA/NA»NA» 'FPE ') 

N2 = NAINA 

CALL READMA(FTT(IST+N2bNAiNA»li .FALSE. » '(3D20. 13) ') 

CALL URITMA(FTT(IST+N2bNAiNA»b 'GPE ') 

CALL READMA ( FTT ( I ST+N2+NA b b b NA > » TRUE ♦ » ' ( 3D20 . 1 3 ) ') 

CALL HR I TMA ( FTT ( I ST+N2+NA b b b NA i ' HPE ') 

CALL READMA ( FTT ( 1ST FN2+NAI2 b b b b . TRUE . » ' ( 3D20 ♦ 1 3 ) ') 

CALL URITMA(FTT<IST+N2+2INAbbl»b'JPE ') 

1ST = 1ST +N2+3INA+2 
C 

100 IFi.NOT. LX(27)) GO TO 110 
URITE(3»9030> 

9030 FORMAT (/» ' RUDDER CONTROL PILOT MODEL MATRICES') 

READ(lrlGl) NA 
IFPR(l) ^ NA 
IFFR(2) = NA 
1FPR<3) = 1ST 

CALL READMA (FTT (I ST) »NA»NA»NA».TRUE.» '(3D20.13) ') 

CALL WRITMA(FTT(I3T)»NA»NA/NA/'FPR ') 

N2 = NAINA 

CALL RE ADMA ( FT T ( I ST+N2 b NA / N A * b » F ALSE . / ' £ 3D20 . 1 3 ) ') 

CALL HRITMA(FTT(ISTFN2bNA/NAfb'GPR ’) 

CALL READMA OFTT < IST+M2+NA) »1 f 1 »NA» >TRUE* » ' (3D20. 13) ') 

CALL WRITMA(FTT(IST+N2FNAbbbNA»'HPR ') 

CALL READMA (FTT ( 1ST +N2+2INA) r 1» b li .TRUE. » ' (3D20. 13) ') 

CALL URITMA ( FTT < 1ST FN2+2INA b b b lr ' JPR ') 

1ST = 1ST + N2+3INA+2 
C 

110 IFi.NOT. LX (28)) GO TO 120 
WRITE(3»9040) 

9040 FORMAT (/'»' THRUST CONTROL PILOT MODEL MATRICES') 

READ(hlOl) NA 
IFFT(l) * na 
IFPT (2) * NA 
IFPT (3) = 1ST 

CALL READMA(FTT(ISTbNA/NArNA».TRUE.»'(3D20.13) ') 

CALL HRITMA(FTT(IST)iNAfNAiNA»'FPT ') 

N2 » NAINA 

CALL READMA ( FTT ( I STFN2 b NA » NA 1 1 > .FALSE . » ' (3020.13) ') 

CALL WRITHA(FTT(IST+N2bNA»NA»li'GPT ') 

CALL RE ADM A ( FTT < I ST+N2+NA b b b NA / . TRUE . i ' ( 3D20 . 1 3 ) ’) 

CALL UR1 TMA ( FTT ( 1ST FN2+NA) » b bNA» ' HPT ') 

CALL READMAfFTT ( 1ST +N2+2INA ) i b li b • TRUE * » ' ( 3D20 .13) ') 

CALL NRITMA(FTT ( 1ST +N2+2IMA ) » b li b ' JPT ') 

1ST = IST+N2+NAI3+2 
C 


120 CONTINUE 
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IWK<3> = 1ST 
C 

C FTTC INPUT PARAMETERS 

C 

READ ( 1 » 101 > NSAMPZ.NSAHPP 
URITE(3,108) NSAMPZ » NS AMPP 

108 FORMAT < /» ' SAMPLING MULTIPLES OF INTEGRATION TINE FOR THE'./ 

COMPENSATOR AND PILOT'. 10X.2I5) 

URITE<3»109) 

10? FORMAT ( /, ' AIRCRAFT TRIM VALUES') 

CALL READHA(US! 11 ) . 1 . 1 »NZF4. .TRUE. » ' (5E10.3) ') 

CALL WRITMA(USdl).l.l.NZ.'YTRM') 

CALL WRITHA (US( 1 1+MZ) » 1 .1 .4. 'UTRM' ) 

C 

WRITE( 3 . 107 ) 

107 FORMAT!/, SENSOR NOISE STANDARD DEVIATIONS') 

CALL READMA(US(26),1,1,NZ,.TRUE.,'!5E10.3) ') 

CALL NRLTMA(US(26),1»1,NZ,'RZ ') 

C 

READ! 1.101) NTRAJP 
NRITE!3,106) 

106 FORMAT!/,' TRAJECTORY PARAMETERS') 

CALL READMA(US(41),1,1»NTRAJP».TRUE., '(5E10.3) ') 

CALL NR I TMA ( US ( 4 1 ) . 1 . 1 , NTRAJP . ' TRA J ' ) 

C 

HRITE(3»105) 

105 FORMAT!/,' CONTROL UNIT CONVERSIONS (COMP. TO SIMII)') 

CALL READMA(US(56),1,1,4,.TRUE.,'!5E10.3) ') 

CALL WRITMA(US<56) , 1 » 1 » 4 » 'UNIT') 

C 

C INITIALIZE FOR SENSOR NOISE 

N - 0 

ST = GRAND! N) 

C 

RETURN 

END 
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SUBROUTINE SEHSOR(Z»ZTRIH»ZC»EC»RZfNZ> 

IMPLICIT DOUBLE PRECISION (A-H,0-Z> 

DIMENSION Zd)>ZTRIM(l)>ZCd)rEC(l)fKZd) 

COMMON /'DRVOUT/' F(13)»DF(13> 

COMMON /5IM0UT/ AMCH rQBARr GMA r DEL » UB > VB t HB » VEAS * OCAS 
COMMON /SIMACC/ AX»AYiAZ»ANX»ANY»ANZfAN 
C 

C Z(I) - UXdO+I) 

C 

C ZTRIM(l) - USdO+I) 

C ZCd) = UXC25+I) 

C ECd) ^ UX(40+I> 

C RZd) = US (25+1) 


C 

C 

C OBTAIN SPECIFIC MEASUREMENTS HERE FROM SIMULATION 

C OUTPUTS 

C START SPECIFIC ZAPO CODE HERE (SEE FTTC RPT SEC 5) 

C UNITS-- FT/SEC AND .01 RAD» RAD/SEC. ETC. FOR 

C THE LINEAR ZAPO MODEL 

C 

C Zd) -- VERT. ACC. 

C Z(2) -- ANGLE OF ATTACK 

C Z(3) — PITCH RATE 

C Z<4) -- FITCH ANGLE 

C Z(5) -- MACH 

C 
C 

C FORM SENSOR MEASUREMENTS IN SIHII UNITS 

C 


Zd) ^ ANZ + C0S(F(8))*C0S(F(10>) 
Z(2) = F(6) 

Z(3) = F(3) 

Z(4) - F(8) 

Z(5) = AMCH 


C FORM ERROR SIGNAL IN COMPENSATOR UNITS 

ECd) = (Zd) -ZTRIM(l) - ZC(1))*32.174 
EC(2) = (Z(2) -ZTRIM(2) - ZC(2))*100. 

EC(3) = (Z(3) -ZTRIH<3) - ZC(3))*100. 

C TURN GAMMAC INTO THETAC 

EC(4) = (Z(4) -ZTRIM(4) - (ZC(4)+ZC(2)))*100. 

EC ( 5) = ( Z ( 5) -ZTRIM(5) - ZC(5))I1000. 

C END SPECIFIC ZAPO CODE HERE 

C 

C ADD MEASUREMENT NOISE ( MUST HAVE COMPENSATOR UNITS!) 

DO 10 1=1 »NZ 
ST = GRAND(N) 

C TEMP HOD TO SEE NOISE IN UX(61 )-UX( 65) 


10 Z(I) ^ Zd) + STtRZ(I) 
RETURN 
END 
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SUBROUTINE TRAJGN(TRAJV,ZC) 

IMPLICIT DOUBLE PRECISION (A-H»0-Z) 
DIMENSION TRAJVUbZCU) 


C 

C TRAJGN GENERATES THE TRAJECTORY COMMANDS AND STORES THE 
C COMMANDED OUTPUTS IN ZC 
C 

COMMON /DRVOUT / F<13bDF(13) 

COMMON /SETICS/ FI(13) 

C 

C THE FOLLOWING CODE MUST BE MODIFIED FOR DIFFERENT 

C CLASSES OF MANEUVERS 

C 

C START SPECIFIC ZOOM AND PUSH OVER CODE 

C (SEE FTTC FINAL REPORT SECTION 5,1) 

C 

C TWO CONSTANTS AND THE DESIRED BIAS VALUES ARE THE 

C ? TRAJECTORY VALUES 

C TRAJV(l) = AZO 

C (2) = ALPHAO 

C (3) = QO 

C (A) = GAHHAO 

C (5) = MO 

C (6) » TFINAL 

C (7) - GIR0/CW2 

C 


ZCU) = TRAJV(l) 

ZC(2) = DALPHAIF ( 1 )/'TRAJV<6) +TRAJV(2) 

ZC( 3) = TRAJV(3) 

ZC<4) = -GAHAIC4F ( 1 )/TRAJV(6) + TRAJV(4) E GAMAIC 
ZC(5) = DSQRT(TRAJV(7)*ZC(4)«2 E TRAJV(5)U2) 
RETURN 
C 

C DO TRAJECTORY PARAMETER SETUP 

C 

ENTRY TRAJST 
GAMAIC - FI(8)-FI(6) 

DALFHA = TRAJV(2> - FI (6) 

C 

C STOP ZAPO SPECIFIC CODE 

C 

RETURN 

END 
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